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UPPER  MANTLE  COMPRESSIONAL  WAVE  VELOCITY  BENEATH  IRAN 


by  Keith  Priestley 


Abstract 

The  upper  mantle  compressional  wave  velocity  beneath  the  Iranian 
Plateau  has  been  determined  by  measuring  two  station  interval  velocities 
from  first  arrivals  in  the  distance  range  2°  to  18°  reported  in  the  Bulletin 
of  the  International  Seismological  Centre  (ISC).  These  results  show  that 
the.Pn  velocity  is  near  8.1  km/sec  beneath  the  Zagros  Mountains  and  most 
of  the  northern  part  of  the  Iranian  Plateau  and  is  not  as  variable  as  re¬ 
ported  in  past  studies.  A  lower  velocity  of  about  7.8  km/sec  occur 
beneath  northwestern  Iran.  The  magnitude  of  various  sources  of  error 
are  discussed.  Teleseismic  P-wave  delays  are  used  to  estimate  the  crustal 
delays  and  the  apparent  Pn  velocities  are  corrected.  Additional  P„  velocity 
estimates  are  made  using  digital  recordings  of  well  located  aftershocks  of 
the  1978  Tabas  earthquake  in  northeastern  Iran  from  the  Iranian  Long  Pe¬ 
riod  Array.  Pn  velocity  from  the  aftershock  seismograms  confirm  the  two 
station  Pn  velocity  measurements  from  the  ISC  data.  Using  ISC  travel 
times  from  events  with  focal  depths  constrained  by  waveform  modelling 
proved  to  still  have  too  great  a  scatter  to  refine  the  Pn  velocity  estimates 
over  previous  single  station  measurements  in  this  region. 


Introduction 


This  study  attempts  to  refine  the  measurement  of  upper  mantle  compressional  velocity 
beneath  the  high  topographic  region  of  Iran.  Accurate  Pn  velocity  estimates  are  an  important 
parameter  in  locating  small  magnitude  events  in  this  region  since  most  such  events  are  located 
with  a  few  seismographs  at  regional  distance  ranges.  Several  previous  studies  have  estimated 
the  Pn  velocity  in  this  region  but  the  results  are  conflicting.  For  example  Chen  et  al  (1980), 
Kadinsky-Cade  et  al  (1981)  and  Hearn  &  Ni  (1994)  all  find  low  upper  mantle  velocity  in 
the  region  of  the  Turkish-Iranian  border  NW  of  the  station  TAB  (Fig.  1),  but  SE  of  TAB 
in  Iran  Chen  et  al  (1980)  and  Asudeh  (1982a)  find  the  Pn  velocity  to  be  about  7.8  km/sec 
while  Kadinsky-Cade  et  al  (1981)  find  the  Pn  velocity  to  be  8.6  km/sec. 

Seismologist  use  the  term  Pn  for  two  distinctly  different  seismic  phases.  In  controlled 
source  refraction  experiments  the  term  Pn  is  used  to  denote  the  headwave  propagating  im¬ 
mediately  beneath  the  Moho.  The  apparent  velocity  of  this  phase  in  the  case  of  plane, 
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horizontal  layers  is  equivalent  to  the  upper  mantle  velocity  directly  directly  beneath  the 
Moho.  Depending  on  the  upper  mantle  velocity  gradient  (Hill,  1971),  the  Pn  headwave  can 
be  observed  to  a  distance  of  several  hundred  kilometers,  after  which  it  is  typically  replaced 
by  sets  of  delayed  arrivals  with  a  similar  apparent  velocity  (Him  et  al,  1973).  Earthquake 
seismologists  use  the  term  Pq  to  denote  the  first  arrival  from  earthquakes  in  the  distance 
range  from  about  200  km  to  as  great  as  2000  km.  Beyond  about  2000  km  the  Pj,  phase  is 
replaced  as  the  first  arrival  by  direct  P  turning  below  the  410-discontinuity,  but  Pj,  can  often 
be  observed  as  a  secondary  arrival  to  greater  distance.  This  regional  distance  range  seismic 
phase  has  nearly  constant  phase  velocity  over  this  entire  distance  range.  This  second  sense 
of  the  term  Pn  is  used  in  this  work. 

In  an  axea  such  as  Iran,  the  determination  of  upper  mantle  velocities  must  depend  on 
seismograms  from  earthquakes.  However  it  is  well  known  that  earthquakes,  especially  smaller 
events  (mb  <  5),  are  poorly  located  (Asudeh,  1983).  Many  regions  in  Iran  have  high  seis¬ 
micity  but  the  seismograph  density  is  too  low  to  permit  the  simultaneous  solution  for  the 
velocity  and  hypocentral  parameters.  In  this  study  I  have  examined  three  ways  of  improv¬ 
ing  the  estimates  of  the  upper  mantle  compressional  velocity;  (a)  using  two  station  interval 
velocity  measurements  from  regional  arrival  times  reported  by  the  Bulletin  of  the  Interna¬ 
tional  Seismological  Center  (ISC);  (b)  using  digital  recordings  of  well  located  aftershocks  of 
the  Tabas  earthquake;  and  (c)  using  ISC  arrival  time  data  from  moderate  events  whose  focal 
depths  are  constrained  by  waveform  modelling. 

Interval  Velocity  Estimates  Using  ISC  Arrival  Times 


METHOD 

The  two  station  interval  velocity  method  is  a  modified  form  of  the  time-term  method 
of  Scheidegger  &  Willmore  (1957)  and  Willmore  &  Bancroft  (1960).  Several  studies  have 
used  this  method  to  estimate  uppermost  mantle  velocity  from  arrival  times  reported  by  the 
ISC  (Asudeh,  1982a;  Beghoul  &  Barazangi,  1989)  or  other  networks  (Haines,  1979,  1980; 
Kayal  &  Smith,  1984;  Kayal  &:  Reena  De,  1987).  Asudeh  (1982a)  used  ISC  arrival  time 
data  from  Iranian  seismographs  to  determine  Pn  velocity  along  several  paths  to  a  precision 
of  considerably  better  than  0.1  km/sec. 

If  two  seismographs  Si  and  S2  lie  on  the  same  great  circle  path  as  an  earthquake  (Fig. 
2),  then  the  apparent  upper  mantle  velocity,  c,  between  the  two  stations  is  given  by 


where  Ai^2  is  the  epicentral  distance  and  Ti^2  is  the  arrival  time  at  the  seismographs  1  and  2, 
respectively.  If  A  is  in  the  range  of  a  few  degrees  the  curvature  of  the  Earth  may  be  ignored 
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and  if  the  crustal  delays  are  the  same  beneath  both  stations,  the  uppermost  mantle  velocity, 
V,  between  the  two  stations  is  equal  to  the  apparent  velocity. 


where  A  and  T  are  the  difference  in  epicentral  distance  and  arrival  time,  respectively.  Errors 
in  the  focal  depth  and  origin  time  are  eliminated  by  taking  the  arrival  time  difference  in 
(1).  If  the  crustal  thickness  or  velocity  varies  beneath  the  two  seismographs  then  different 
apparent  velocities  will  be  observed  along  the  forward  and  reversed  path.  In  this  case  the 
true  upper  mantle  velocity  will  be  approximately  the  average  of  the  forward  and  reversed 
apparent  velocities. 

Earthquakes  in  Iran  are  known  to  be  poorly  located  by  the  ISC  (Asudeh,  1983).  Location 
errors  will  map  into  errors  in  the  measured  velocity.  These  probably  account  for  much  of 
the  scatter  in  travel  time  curves  derived  from  a  single  station  when  the  ISC  locations  axe 
used,  even  when  the  arrival  times  are  determined  from  the  original  records  (Chen  et  al, 
1980;  Kadinsky— Cade  et  al,  1981).  Asudeh  (1982a)  and  Beghoul  &  Barazangi  (1989)  have 
shown  that  the  two  station  method  minimized  the  effect  of  location  errors.  Following  Asudeh 
(1982a),  if  E'  (Fig.  2)  is  the  epicenter  reported  by  the  ISC  and  E  is  the  true  epicenter,  then 
EE'  (with  components  6x,  6y,  and  Sz)  is  the  mislocation  vector.  The  effect  of  6x  and  Sz 
cancel  in  taking  the  arrival  time  difference  T2  —  Ti  and  distance  difference  A2  —  Ai,  and  only 
the  component  Sy  perpendicular  to  the  great  circle  path  will  affect  the  velocity  calculation. 
From  Figure  2 


Sy 

sm  ®  ^  •  (3) 

Small  errors  in  the  location  and  interstation  travel  time  leads  to  errors  in  the  velocity,  Sv. 
From  (2)  and  (3) 

A  cos  Tf 

^  ~  ,  (4) 


thus 


Sv 


—A  cos  '1? 
2^2 


•  ST  - 


A 

T 


sin  ^  •  SQ 


(5) 


For  magnitude  4.5  to  5  events  in  Iran  the  mislocation  vector  EE'  is  likely  to  be  10  to  20  km 
(Asudeh,  1983).  The  crust  beneath  the  Iranian  plateau  is  40  to  50  km  thick  (Asudeh,  1982b) 
therefore  Aj  is  about  200  km.  For  Ai  greater  than  about  400  km,  which  is  more  typical  of 
the  events  used  in  this  study,  the  error  introduced  by  the  second  term  on  the  right-hand 
side  of  (5)  is  less  than  0.01  km/s  (Asudeh,  1982a)  and  can  be  ignored  and  (5)  reduces  to 
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(6) 


6v 


—A  cos  0 
jJi 


ST  . 


Thus  errors  in  the  epicenter  have  a  small  effect  on  the  velocity  calculation;  errors  in  the  focal 
depth  and  origin  time  cancel  entirely.  In  the  absence  of  clock  errors,  ST  reflects  timing  errors 
which  should  be  normally  distributed  and  hence  should  be  reduced  by  repeated  measurement. 
Also,  the  magnitude  of  Sv  is  inversely  proportional  to  T,  so  large  A  which  give  large  T  lead 
to  more  accurate  velocity  estimates.  For  A  700  km  and  Pj,  «  8.1  km/sec  which  is  typical 
in  this  study,  a  one-second  reading  error  gives  less  than  0.1  km/sec  error  in  the  Pn  velocity. 

DATA 


Data  for  the  two  station  interval  velocity  measurements  have  been  taken  from  the  ISC 
CD-rom.  This  includes  data  from  January  1,  1964  to  August  31,  1987.  Arrival  time  data 
for  the  sixteen  stations  in  Iran,  eastern  Iraq,  and  southern  Azerbijan  shown  in  Figure  1  were 
used  (Table  1).  All  crustal  events  for  which  (a)  the  ISC  had  used  at  least  25  stations  in  the 
location,  (b)  the  stations  lay  within  the  epicentral  distance  range  2°  to  18°,  and  (c)  whose 
travel  time  residual  was  less  than  5  seconds  were  chosen.  The  first  restriction  ensured  that 
the  epicenters  were  reasonably  well  located.  In  most  cases  the  number  of  locating  stations 
was  much  greater  than  25.  The  second  restriction  ensured  that  the  first  axrival  was  in  the 
Pn  distance  range.  Timing  errors  are  common  in  the  ISC  arrival  data.  Many  of  these  are 
obvious  10  and  60  second  notation  errors.  The  third  restriction  removes  the  affect  of  these 
large  timing  errors.  To  obtain  a  large  number  of  arrival  time  data  the  requirement  that  the 
event  and  two  stations  lie  on  the  same  great  circle  path  was  relaxed.  Both  Asudeh  (1982a) 
and  Beghoul  &  Barazangi  (1989)  allowed  the  opening  angle  0  (Fig.  2)  to  be  6°.  In  this  study 
0  equal  7°  was  used  in  selecting  events.  The  increased  error  from  using  this  larger  value 
of  0  is  probably  insignificant  compared  to  that  from  using  6°  as  in  Beghoul  &:  Barazangi 
(1989)  and  0  =  7°  seems  to  correspond  to  what  Asudeh  (1982a)  actually  used  to  obtain  the 
events  he  shows  in  his  Table  1.  The  selection  criteria  will  be  discussed  further  below  when 
the  effects  of  errors  are  assessed. 

RESULTS 

The  routines  for  selecting  the  data  from  the  ISC  file  and  for  computing  the  two  station 
interval  velocity  were  first  tested  to  reproduce  the  results  in  Table  1  of  Asudeh  (1982a). 
Other  than  specifying  a  maximum  opening  angle  of  6°,  Asudeh  (1982a)  does  not  cite  his 
criteria  for  chosing  his  data.  Using  the  portion  of  the  ISC  data  between  July  26,  1967  and 
October  2,  1976,  a  number  of  test  selections  were  made  to  determine  the  criteria  Asudeh 
used.  To  select  the  seven  events  in  his  Table  1  required  0  =  7°,  travel  time  residuals  less  than 
5.5  seconds,  and  ISC  focal  depths  to  57  km.  Many  ISC  focal  depths  for  events  in  Iran  are 
sub-Moho  but  there  is  no  reliable  evidence  for  mantle  events  except  in  the  Makran  (Jackson 
&  Fitch,  1981).  There  are  a  large  number  of  events  which  occur  during  this  time  period 
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which  also  satisfy  these  selection  criteria  but  which  Asudeh  did  not  use  and  he  does  not 
comment  why  the  events  in  his  Table  1  were  chosen  and  the  other  events  rejected.  Using  the 
arrival  time  data  from  the  CD-rom,  phase  velocity  values  were  determined  to  within  ±0.02 
km/sec  for  six  of  the  seven  values  shown  in  Table  1  of  Asudeh  (1982a),  but  the  apparent 
velocity  of  the  October  2,  1976  event  was  nearly  0.4  km/sec.  The  arrival  time  difference  for 
this  event  in  Table  1  of  Asudeh  (1982a)  is  wrong  by  4  seconds. 

The  two  station  interval  velocity  measurements  are  show  as  a  function  of  position  and 
direction  in  Figure  3.  These  values  plus  their  average  error  are  tabulated  in  Table  1.  North¬ 
ern  Iran  is  well  sampled,  however,  the  values  for  central  Iran  all  depend  on  the  WWSSN 
station  SHI.  The  apparent  Pj,  velocity  in  central  and  NE  Iran  is  8.1±0.1  km/sec.  Beneath 
the  T^esh  and  Albroz  Mountains  along  the  Iranian  shore  of  the  Caspian  Sea  and  in  NW 
Iran  near  the  Turkish  border,  the  Pj,  velocity  is  less  than  8.0  km/sec.  The  apparent  P^ 
velocity  measurements  in  Figure  3  suggest  less  variability  in  the  upper  mantle  compressional 
wave  velocity  than  found  in  previous  studies.  The  velocity  average  errors  in  Table  2  are 
significantly  larger  than  those  reported  in  Asudeh  (1982a).  Ignoring  the  error  in  his  Table 
1,  the  lower  errors  he  reports  seem  to  arise  from  the  way  he  has  culled  the  data. 

ACCESSMENT  of  ERRORS 

The  effect  of  epicentral  mislocations  has  been  discussed  above.  However,  there  are  several 
other  potential  sources  of  error  in  the  procedure  outlined  above. 

Effect  of  Spherical  Wave  Front  For  0  greater  than  zero,  the  true  curved  wavefront  arrives 
at  the  distant  station  at  a  later  time  than  if  the  wavefront  was  planar  as  assumed.  Beghoul 
&  Barazangi  (1989)  show  that  this  difference  in  distance  is  given  by 

.  HH'  =  A2  -  Ai  sin^  0  -  cos  ©  -  (Ai)2  sin"*©  (7) 

and  that  for  0  equal  6°  and  an  interstation  distance  of  14°  HH'  is  7.4  km  which  leads  to 
less  than  a  0.5%  error  in  the  estimated  P„  velocity.  Thus  neglecting  the  wavefront  curvature 
leads  to  insignificant  error. 

Effect  of  Variation  in  0:  There  is  a  trade-off  in  choosing  0.  To  obtain  a  large  number 
of  events  to  make  the  interval  velocity  estimates,  larger  0  is  desirable.  However,  larger  0 
couples  with  epicentral  mislocation  to  produce  larger  errors  in  the  apparent  P^  velocity. 
Beghoul  &  Barazangi  (1989)  constrained  0  to  6°  and  from  tests  of  the  data  in  Table  1  of 
Asudeh  (1982)  he  appears  to  have  let  0  be  at  least  7°.  Figure  4  shows  the  combined  effect  of 
varying  0  and  epicentral  mislocation.  The  magnitude  of  the  error  decreases  with  increasing 
distance.  Figure  4a  shows  the  error  for  0  of  14°  and  an  epicentral  mislocation  of  28  km 
(20  km  perpendicular  to  the  interstation  line  and  20  km  in  line  with  the  interstation  line). 
Figure  4b  assumes  the  same  epicentral  mislocation,  but  is  for  0  of  7°  and  is  the  data  used 
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in  measuring  the  apparent  velocities  shown  in  Figure  3.  For  this  0  and  for  an  assumed 
epicentral  mislocation  of  28  km,  which  should  be  representative  for  the  Middle  East  events 
used  in  this  study,  the  error  is  less  than  1.5%. 

Effect  of  Lid  Velocity  Gradient:  The  propagation  of  the  mantle  lid  P-wave  is  primarily  con¬ 
trolled  by  the  lid  compressional  wave  velocity  gradient  including  the  pseudo-velocity  gradient 
caused  by  Earth  sphericity.  Braile  &  Smith  (1975)  used  reflectivity  synthetic  seismograms 
to  show  the  effect  of  the  velocity  gradient  directly  below  the  Moho  on  Pn  propagation.  Even 
in  the  presence  small  negative  velocity  gradient  in  the  lid,  the  pesudo-gradient  as  a  result 
of  curvature  results  in  turning  rays  in  the  lid  and  Pn  propagates  to  distances  of  400  to  700 
km.  Pn  is  not  observed  only  in  cases  where  the  negative  gradient  is  enough  to  off-set  that 
resulting  from  curvature.  Menke  &  Richards  (1980)  modelled  teleseismic  Pn  as  a  sum  of 
whispering  gallery  waves  which  propagate  in  a  waveguide  composed  of  the  high  velocity 
mantle  lid  overlying  the  upper  mantle  low  velocity  zone.  The  velocity  gradient  of  the  lid, 
including  the  pseudo-gradient  resulting  from  Earth  curvature,  is  the  most  important  pa¬ 
rameter  controlling  propagation  and  both  crustal  structure  and  upper  mantle  velocity  are 
of  secondary  importance.  Lid  velocity  structures  which  have  weak  velocity  gradients  have 
nearly  constant  apparent  velocity  with  increasing  range.  The  observation  that  Pn  falls  along 
a  linear  travel-time  branch  suggests  that  velocity  gradient  in  the  mantle  lid  is  low  and  hence 
the  Pn  velocity  found  at  larger  distances  is  representative  of  the  true  upper  mantle  material 
speed. 

The  two  station  method  assumes  that  both  stations  record  the  same  phase.  If  a  strong 
positive  velocity  gradient  exist  in  the  mantle  lid  then  the  two  station  Pn  velocity  estimate 
would  depend  on  the  epicentral  distance  of  the  two  stations  since  the  near  station  would 
be  recording  a  shallow  turning  phase  while  the  far  station  recorded  a  deeper  turning  phase. 
To  estimate  the  strength  of  the  lid  gradient  beneath  the  Middle  East,  travel  times  for  all 
stations  in  the  distance  range  2°-25°  were  plotted  as  a  reduced  travel  time  curve.  Travel 
time  curves  for  four  stations  axe  shown  in  Figure  5.  These,  and  those  for  the  other  twelve 
stations  not  shown,  indicate  that  the  first  arrival  in  the  distance  range  from  about  2°  to 
about  18°  fall  on  a  linear  travel  time  branch.  Based  on  this  observation,  only  events  more 
than  200  km  from  the  near  station  and  less  than  1800  km  from  the  far  station  were  used  in 
computing  the  apparent  Pn  velocity. 

Table  3  gives  the  intercept  times  for  the  ISC  arrival  time  data  for  each  of  the  stations. 
In  computing  these  intercept  times,  only  arrivals  from  events  in  the  distance  range  250  to 
1600  km  with  residuals  less  than  5  seconds  were  used. 

Effect  of  crustal  delay  times  beneath  the  stations:  The  source  crustal  delay  term  is  common 
to  the  arrival  times  at  both  stations  and  therefore  cancels  when  taking  the  time  difference  in 
the  interval  velocity  calculation.  In  the  case  of  a  crust  of  uniform  thickness  and  velocity,  the 
apparent  P„  velocity  from  the  two  station  interval  velocity  measurement  would  be  equivalent 
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to  the  upper  mantle  compressional  wave  velocity.  In  the  case  where  the  crustal  structure 
changes  between  receivers  but  where  reversed  apparent  velocity  estimates  are  available,  the 
upper  mantle  velocity  is  approximately  the  average  of  the  two  way  measurements.  However, 
variations  in  crustal  thickness  or  crustal  velocity  in  the  vicinity  of  the  receivers,  neither  of 
which  are  known  for  Iran,  can  combine  to  produce  time  delays  which  would  bias  the  estimate 
of  the  Pn  velocity.  Even  though  the  crustal  structure  of  Iran  is  poorly  known  we  can  make 
some  estimate  of  the  magnitude  of  the  effects  of  the  crustal  delays.  The  magnitude  of  the 
effect  of  crustal  structure  variation  depends  on  the  path  length  between  the  two  stations.  For 
a  path  length  of  250  km,  replacing  5  km  of  crust  of  velocity  6.1  km/sec  with  upper  mantle 
material  of  velocity  8.1  km/sec  changes  the  estimated  Pn  velocity  from  the  true  velocity  of 
8.10  km/sec  to  8.18  km/sec  but  for  a  path  length  of  750  km  a  5  km  thinning  of  the  crust 
only  changes  the  estimated  P^  velocity  to  8.12  km/sec. 

One  estimate  of  the  variation  in  crustal  delay  terms  exist  through  comparing  relative 
teleseismic  delay  times  across  Iran.  Figure  6  (Table  4)  is  a  plot  of  the  teleseismic  delay  times 
at  15  seismographs  relative  to  the  station  SHI  in  the  Zagros  Mountains.  These  residuals 
were  determined  by  finding  the  mean  of  the  difference  in  the  P-wave  residual  at  a  recording 
station  minus  the  residual  at  station  SHI.  Station  SHI  was  chosen  because  it  had  the  largest 
number  of  recordings.  Only  events  greater  than  35°  distance  and  which  were  not  emergent 
arrivals  were  used.  With  one  exception.  Figure  6  shows  a  very  consistent  pattern.  Except  for 
the  station  MSH,  all  the  stations  in  eastern  Iran  show  relative  residuals  of  about  0.55±0.15 
second.  Four  stations  MSH,  MHI,  MAIO,  and  MUI  are  located  close  together  at  Mashhad. 
MHI  and  MAIO  are  co-located  and  have  similar  relative  travel  time  residuals.  MUI  and 
MSH  are  located  about  15  km  from  MHI  and  MAIO,  and  about  2  km  from  each  other. 
MUI  has  a  relative  residual  similar  to  that  at  MHI/MAIO  but  the  delay  at  MSH  is  much 
greater.  Dziewonski  &  Anderson  (1983)  examined  teleseismic  P-wave  residuals  from  the 
ISC  arrival  time  data  and  found  that  MSH  showed  much  larger  delays  than  nearby  stations. 
MSH  also  produces  anomalously  low  velocity  estimates  from  the  first  arrival  of  regional 
earthquakes.  Without  knowing  more  details  of  the  site  condition  at  MSH,  it  is  impossible  to 
determine  whether  the  station  has  a  large  relative  residual  compared  to  nearby  stations  due 
to  site  condition,  such  as  being  located  above  a  sedimentary  basin,  or  whether  other  factors 
contribute  to  the  larger  delay  at  MSH.  Since  MSH  produced  no  interval  velocity  estimates 
which  fell  in  the  Pn  velocity  range  and  since  MSH  shows  more  than  a  half  second  larger 
delay  compared  to  the  nearby  stations,  data  from  MSH  has  not  been  used  in  this  study. 
There  is  about  a  half  second  relative  delay  from  the  northeastern  Iranian  stations  compared 
to  the  Zagros  Mountains  stations.  If  this  is  solely  do  to  crustal  thinning  this  would  imply 
that  the  crust  beneath  the  Zagros  is  about  5  km  thinner  than  the  crust  beneath  northeast 
Iran.  However,  both  lower  crust  or  upper  mantle  velocity  beneath  NE  Iran  could  contribute 
to  their  delays. 

The  three  stations  in  northern  Iran  south  of  the  Caspian  Sea  (TEH,  ILPA,  SRI)  show 
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similar  teleseismic  delays  of  about  1  second  relative  to  SHI.  If  the  crustal  velocities  in  the 
two  regions  are  similar,  the  crust  beneath  the  Albroz  Mountains  is  about  10  km  thicker  than 
the  crust  beneath  the  Zagros  Mountains.  The  station  TAB  in  northwest  Iran  and  the  nearby 
station  LNK  in  southeastern  Azerbijan  show  similar  delays  of  about  1.45  seconds,  although 
there  are  few  arrival  times  for  the  LNK  station.  TAB  is  in  a  region  of  recent  volcanism  and 
this  has  previously  been  cited  as  a  possible  explanation  for  the  low  P„  velocities  reported  in 
this  area  (Chen  et  al,  1980;  Kadinsky-Cade  et  al,  1981;  Asudeh,  1982a).  This  might  also 
explain  part  of  the  large  teleseismic  delay  observed  at  TAB  relative  to  SHI.  However  the 
delay  at  TAB  is  identical  to  that  observed  at  LNK  200  km  ENE. 

LNK  is  the  only  station  in  this  region  where  the  crustal  structure  beneath  the  station  is 
well  known.  Mangino  k  Priestley  (1996)  have  analyzed  broadband  teleseimic  waveforms  from 
a  digital  seismograph  co-located  with  the  analogue  seismograph  which  reports  arrival  times 
to  the  ISC.  To  determine  the  crustal  structure  they  have  inverted  radial  receiver  function 
data  using  the  routine  of  Ammon  et  al  (1990),  then  forward  modelled  the  radial  receiver 
function  to  assess  which  features  of  the  inversion  model  were  significant.  The  results  of 
this  analysis  for  LNK  are  summarized  in  Figure  7.  The  inversion  results  (Fig.  7a)  imply 
the  presence  of  a  low  velocity  (Vp  <  4.8  km/s)  upper  most  crust  (0-6  km),  a  large  (>2.0 
km/s)  velocity  contrast  between  6-12  km  depth,  a  high  velocity  (Vp  7.0  km/s)  mid-  to  lower- 
crust  and  a  complex  crust-mantle  boundary.  Between  8- 12s  the  amplitude  of  the  tangential 
component  is  about  equal  to  that  of  the  radial  and  there  is  a  negative  gradient  between 
20-25  km  depth  which  is  producing  the  synthetic  radial  motion  at  this  time.  The  forward 
modelling  result  (Fig.  7b)  shows  that  the  the  main  features  of  the  radial  receiver  function 
can  be  fit  with  a  much  more  simple  crustal  velocity  model.  Thus  the  LNK  teleseismic  delay 
is  largely  the  result  of  the  low  velocity  upper  crust.  This  is  not  likely  to  be  representative  of 
the  structures  in  the  Zagros  Mountains  and  the  Iranian  plateau. 

The  actual  crustal  thickness  is  not  important.  If  the  teleseismic  delays  in  Figure  6  are 
the  result  of  varying  crustal  thickness  or  crustal  velocity  and  not  from  varying  upper  mantle 
velocity,  then  a  crude  correction  of  the  apparent  Pn  velocities  can  be  made  by  converting 
the  teleseismic  crustal  delay  into  an  equivalent  crustal  delay  for  Pn.  Waveform  modelling  of 
regional  seismic  waveforms  of  central  and  southern  Iran  earthquakes  digitally  recorded  on 
stations  of  the  Caspian  Seismic  Network  (Priestley,  work  in  progress),  and  receiver  function 
analysis  for  a  seismograph  station  within  1  km  of  the  Iran  border  in  Turkmenistan  (Mangino 
k  Priestley,  1996)  show  that  the  average  crustal  velocity  in  the  Iranian  plateau  is  about 
6.1  km/sec  which  is  low  for  continental  crust.  Using  this  and  an  average  P„  velocity  of  8.1 
km/sec  gives  a  critical  angle  for  P^  of  48.9°.  The  ratio  of  the  teleseismic  delay  to  the  P^ 
delay  is  inversely  proportional  to  the  ratio  of  the  cosine  of  the  teleseismic  emergent  angle  to 
the  cosine  of  the  P^  critical  angle,  which  in  this  case  is  1.42.  Using  this  factor,  the  teleseismic 
delays  were  converted  to  equivalent  Pn  crustal  delay  and  the  Pn  velocities  shown  in  Figure 
3  were  corrected. 


8 


Figure  8  plots  these  corrected  apparent  velocities  and  the  average  of  the  two-way 
measurements.  Over  most  of  the  Zagros  Mountains  and  the  Iranian  plateau  the  velocity 
is  8.05±0.10  km/sec.  In  northwest  Iran  the  Pn  velocity  is  7.79±0.10  km/sec,  significantly 
lower  then  in  central  Iran. 

Pn  Velocity  Estimates  from  ILPA  Recordings  of  the  Tabas  Aftershocks 

The  September  16,  1978  Tabas— e—Golshan  (Ms  7.7)  earthquake  in  central  Iran  was  fol¬ 
lowed  by  a  large  aftershocks  sequence.  The  University  of  Cambridge  deployed  a  local  seismic 
array  of  nine  stations  in  the  epicentral  region  on  September  28  and  maintained  this  network 
until  October  28.  1560  aftershocks  were  located  and  of  these,  329  are  thought  to  have  un¬ 
certainties  in  the  h3q)ocenter  of  about  1  km  in  the  epicenter  and  2  km  in  focal  depth.  Many 
of  these  have  well  constrained  focal  mechanisms.  The  ILPA  array  («  650  km  distance)  was 
in  partial  operation  for  most  of  this  period  and  ILPA  recordings  of  these  aftershocks  provide 
support  for  the  P^  velocity  from  the  two  station  ISC  data.  Figure  9  is  a  composite  record 
section  consisting  of  seismograms  of  several  Tabas  aftershocks  recorded  across  the  IPLA  ar¬ 
ray.  The  aftershocks  range  in  focal  depth  from  10  to  18  km,  but  have  been  normalized  to  12 
km  depth  (the  median  focal  depth)  using  a  average  crustal  P-wave  velocity  of  6.1  km/sec 
and  a  P^  velocity  of  8.1  km/sec.  The  travel  times  have  been  reduced  by  A/8.1  seconds  and 
the  line  is  for  a  velocity  of  8.1  km/sec.  This  combination  of  sources  and  receivers  averages 
the  Pn  velocity  both  beneath  the  region  of  the  Tabas  sequence  and  the  ILPA  array.  The  8.1 
km/sec  line  closely  matches  the  first  arrivals.  Other  plots  of  multiple  aftershocks  recorded 
at  a  single  ILPA  array  element  suggest  that  the  apparent  Pn  velocity  beneath  the  aftershock 
sequence  may  be  slightly  less  and  the  apparent  P„  velocity  across  the  array  from  a  single 
aftershock  suggest  that  the  Pn  velocity  beneath  the  array  may  be  slightly  higher.  However, 
these  differences  are  probably  not  significant. 

Pn  Velocity  Estimates  from  ISC  Arrival  Times  of 
Large  Well  Recorded  Regional  Events 

A  third  method  of  refining  the  Pn  velocity  beneath  Iran  was  explored  but  proved  unfruit¬ 
ful.  There  have  been  several  studies  of  the  mechanisms  of  moderate  to  large  earthquakes 
in  the  region  of  Iran  using  waveform  modelling  (Baker  et  al,  1993,  Baker,  1993,  Priestley 
et  al,  1994).  Besides  the  fault  plane  solution,  the  initial  P-wave  shape  is  controlled  by  the 
interference  of  direct  P  and  the  depth  phases  pP  and  sP.  The  focal  depths  from  these  studies 
were  used  with  the  ISC  locations  and  the  ISC  first  arrival  times  at  regional  seismograph 
stations  to  derive  single  station  travel  time  curves  corrected  for  source  depth.  Even  with 
improved  focal  depth  estimates,  these  travel  time  curves  showed  too  much  error  to  prove 
useful. 
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Table  1  -  Station  Locations 


Location 

Designation 

Latitude 

(°N) 

Longitude 

(°E) 

Elevation 

(meters) 

Baghdad,  Iraq 

BHD 

33.2744 

44.3558 

32 

ILPA  (site  1),  Iran 

IRl 

35.4162 

50.6887 

1347 

ILPA  (site  2),  Iran 

IR2 

35.6628 

50.8976 

1172 

ILPA  (site  3),  Iran 

IR3 

35.4761 

51.0238 

1106 

ILPA  (site  4),  Iran 

IR4 

35.2387 

50.9012 

1373 

ILPA  (site  7),  Iran 

IR7 

35.7028 

50.6089 

1305 

Kermanshah,  Iram 

KER 

34.3522 

47.1058 

1310 

Kakhk,  Irna 

KHI 

34.1433 

58.6417 

1600 

-  Lenkoran,  Azerbaijan 

LNK 

38.7100 

48.7788 

-2 

Mashhaxd,  Iran  (SR)) 

MAIO 

36.3000 

59.4945 

1150 

Mashhad,  Iran  (WWSSN) 

MHI 

36.3083 

59.4717 

1150 

Mashhad,  Iran  (WWSSN) 

MSH 

36.3111 

59.5878 

987 

Mosul, Iraq 

MSL 

36.3817 

43.1482 

242 

Mashhad  University,  Iran 

MUI 

36.3117 

59.6050 

1000 

Shahrud,  Iran 

SHD 

36.4333 

54.9417 

1500 

Shiraz,  Iran  (WWSSN) 

SHI 

29.6444 

52.5261 

1595 

Safid  Rud,  Iran 

SRI 

36.7583 

49.3833 

243 

Tabriz,  Iran  (WWSSN) 

TAB 

38.0675 

46.3267 

1430 

Tehran,  Iran 

TEH 

35.7378 

51.3856 

1360 

Taghi  Ghanbar,  Iran 

TGI 

32.9611 

59.1933 

1800 
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Table  2  -  Apparent  Velocities 


Path 

Apparent 

No. 

Average 

Path 

Apparent 

No. 

Average 

Velocity 

Error 

Velocity 

Error 

SHD  to  IRl 

8.05 

7 

0.08 

MHI  to  SHD 

7.96 

4 

0.06 

IRl  to  TAB 

7.83 

2 

0.04 

MHI  to  SHI 

8.02 

10 

0.04 

KHI  to  IR2 

8.15 

5 

0.07 

TAB  to  MSL 

8.01 

5 

0.13 

TAB  to  IR2 

7.72 

2 

0.19 

SHD  to  TEH 

7.94 

10 

0.14 

IR2  to  TAB 

7.98 

3 

0.44 

-  - . 

SHD  to  TGI 

8.06 

6 

0.11 

SHD  to  IR7 

8.17 

3 

0.33 

TGI  to  SHI 

8.16 

14 

0.18 

TAB  to  IR7 

8.08 

5 

0.40 

TEH  to  SHI 

8.08 

7 

0.08 

TGI  to  IR7 

8.08 

4 

0.23 

SHI  to  SRI 

7.94 

2 

0.16 

KER  to  SHI 

8.21 

19 

0.18 

SRI  to  SHI 

8.10 

15 

0.08 

SHI  to  KER 

8.16 

27 

0.14 

SRI  to  TAB 

8.09 

4 

0.25 

TAB  to  KER 

7.83 

8 

0.13 

TAB  to  SRI 

8.10 

7 

0.22 

TEH  to  KER 

8.13 

22 

0.16 

SRI  to  TEH 

7.88 

5 

0.20 

MHI  to  KHI 

8.00 

3 

0.06 

TAB  to  TEH 

7.81 

5 

0.13 

KHI  to  MHI 

7.83 

36 

0.12 

TEH  to  TAB 

7.86 

10 

0.14 

KHI  to  SHD 

7.98 

5 

0.11 

TGI  to  TEH 

7.97 

10 

0.12 

SHD  to  KHI 

8.16 

10 

0.09 

TEH  to  TGI 

8.20 

9 

0.16 

KHI  to  SHI 

8.03 

8 

0.09 

MSL  to  BHD 

8.14 

2 

0.25 

KHI  to  TEH 

8.04 

2 

0.01 

KER  to  BHD 

7.99 

3 

0.29 

TEH  to  KHI 

8.24 

2 

0.18 

MSL  to  KER 

8.09 

3 

0.07 

KHI  to  TGI 

8.22 

6 

0.10 

KER  to  MSL 

8.03 

2 

0.05 

TGI  to  KHI 

7.93 

4 

0.19 
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Table  3  -  Pn  Travel  Time  Curve  Intercept  Table  4  -  Teleseismic  Delay  Terms  Relative  to  SHI 
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Figure  Captions 

Figure  1.  Station  locations  shown  on  a  topography  map  of  the  Iran  region. 

Figure  2.  Schematic  representation  of  the  two— station  method,  (a)  Cross  section  showing 
the  P n  ray  paths  to  the  two  seismographs.  The  hase  of  the  crust  is  hatched,  and  the 
vertical  scale  is  exaggerated,  (b)  Map  view  showing  a  two  station  pair  and  defining  the 
opening  angle  0.  Only  events  at  an  epicentral  distance  greater  than  250  km  from  the 
near  station  and  less  than  2000  km  from  the  far  station  are  used,  (c)  An  exaggerated 
schematic  map  view  showing  the  effect  of  the  mislocation  vector  EE'. 

Figure  3.  Apparent  velocity  map  for  the  Iran  region.  The  solid  triangles  denote  the  station 
locations  and  the  arrows  denote  the  direction  of  the  apparent  P„  velocity  estimate.  The 
dot  at  the  tail  of  the  arrow  denotes  the  midpoint  of  the  two  station  path  and  the  value 
beside  the  arrow  is  the  apparent  velocity. 

Figure  4.  Relative  errors  in  the  P^  velocity  versus  average  epicentral  distance  for  two  opening 
angle  (a)  0  =  7°  and  (b)  0  =  14°.  A  location  error  of  28  km  (20  km  perpendicular  to 
the  great  circle  path  connecting  the  stations,  and  20  km  in  line  with  the  stations)  has 
been  assumed.  The  data  shown  in  (a)  are  those  used  to  compute  the  apparent  velocities 
in  Figure  4. 

Figure  5.  Travel  time  curves  for  KER,  SHI,  TAB,  and  KHI.  The  times  have  been  reduced 
by  8  km/sec. 

Figure  6.  Teleseismic  delays  in  seconds  relative  to  the  station  SHI. 

Figure  7.  Receiver  function  structure  for  LNK.  (a)  Inversion  crustal  model,  (b)  Simplified 
crustal  model  from  forward  modelling  of  the  main  features  of  the  receiver  function. 

Figure  8.  Two-way  and  crustal  delay  time  corrected  Pn  velocity  map.  The  one-way,  delay 
time  corrected  Pn  values  are  shown  in  parentheses. 

Figure  9.  Record  section  from  ILPA  recordings  of  the  Tabas  aftershocks.  Depths  of  the  events 
have  been  normalized  to  12  km  and  the  times  have  been  reduced  by  A/8.1  seconds. 
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AVERAGE  LITHOSPHERIC  STRUCTURE  OF  SOUTHERN  AFRICA 

X.  Qiu^,  K.  Priestley,  and  D.  McKenzie 

Bullard  Laboratories,  University  of  Cambridge,  Cambridge  CBS  OEZ,  UK 

SUMMARY 

We  present  a  new  1-D  earth  model  for  the  crust  and  upper  mantle  beneath  the 
Archaean  core  of  southern  Africa.  The  crustal  structure  is  constrained  by  published  seis¬ 
mic  refraction/refiection  data  and  by  modeling  of  teleseismic  receiver  function  data.  The 
mantle  structure  is  constrained  by  travel  time  measurements  for  Pn  and  Sn  waves,  by 
waveform  inversion  of  multi-mode  surface  waves,  and  by  forward  modeling  of  Pn/  wave¬ 
forms  for  regional  earthquakes.  The  42  km  thick  crust  of  southern  Africa  consists  of  four 
layers:  a  2  km  thick  5.05  kms“^  surface  layer,  a  5  km  thick  6.08  kms“^  upper  crust,  a  20 
km  thick  6.30  kms“^  mid-crust,  and  a  15  km  thick  6.73  kms~^  lower  crust.  Below  the 
crust  there  is  an  80  km  thick  upper  mantle  lid.  Vp  and  Vs  are  8.09  and  4.62  kms“^,  re¬ 
spectively,  at  the  top  of  the  lid  and  the  compressional  and  shear  velocity  gradients  through 
the  lid  are  0.0008  and  0.0013  s~^,  respectively.  Below  the  lid  there  is  a  substantial 
shear  wave  LVZ  with  the  shear  wave  velocity  dropping  from  4.72  kms“^  in  the  lower  part 
of  the  lid  to  4.32  km  at  250  km  depth.  The  shear  wave  LVZ  is  required  to  fit  both  the 
Vs„  measurements  and  the  regional  surface  waveforms.  Forward  modeling  of  Pn/  suggests 
that  no  upper  mantle  P-wave  LVZ  exists.  There  is  an  increase  in  the  P-wave  velocity 
gradient  at  125  km  depth  to  0.0015  s“^  and  a  second  increase  at  250  km  depth  to  0.0035 
s“^.  Velocities  and  densities  in  the  seismic  model  of  the  lid  closely  agree  with  velocity  and 
density  estimates  from  geochemical  analysis  of  garnet  peridotite  nodules  from  kimberlites 
on  the  Kaapvaal  Craton,  implying  that  the  estimates  from  kimberlites  are  representative 
of  a  wider  region  of  the  Archaean  core  of  southern  Africa  -  not  just  the  restricted  region 
beneath  where  the  kimberlites  are  found.  A  comparison  of  the  southern  African  seismic 
model  with  seismic  models  for  other  shield  regions  where  nodule  data  are  not  available 
shows  that  the  nodule  results  may  also  be  representative  of  the  upper  mantle  in  those 
regions. 

Key  words:  crustal  structure,  precambrian,  southern  Africa,  upper  mantle 

^Present  address:  South  China  Sea  Institute  of  Oceanology,  Chinese  Academia  of  Sciences,  164  West 
Xingang  Road,  Guangzhou  510301,  P.R.  China 
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1.  INTRODUCTION 


Knowledge  of  the  geological  and  geophysical  characteristics  of  the  continental  lithosphere  and 
upper  mantle  is  a  prerequisite  for  unraveling  continental  evolution.  In  this  study  we  focus  on  the 
seismic  velocity  and  density  structure  of  southern  Africa,  a  region  important  in  understanding 
continental  evolution  for  two  reasons.  First,  xenoliths  from  the  large  number  of  kimberlitic 
extrusions  have  provided  unique  data  on  the  pressure  and  temperature  environment  in  the  upper 
mantle  beneath  southern  Africa  and  estimates  of  upper  mantle  elastic  properties  predicted  from 
the  xenoliths  can  be  tested  with  seismic  measurements.  Second,  because  of  the  vast  amount 
of  geophysical  and  geological  data  accumulated  for  southern  Africa,  this  region  has  become  the 
“Archetype”  of  Archaean  continental  evolution. 

Here  we  analyze  regional  and  teleseismic  distance  range  earthquakes  recorded  on  digital  seis- 
mographs  in  southern  Africa.  We  use  teleseismic  body  waves  in  combination  with  existing  seismic 
refraction  and  reflection  results  to  constrain  the  crustal  velocity  structure,  and  we  use  waveform 
modeling  of  regional  seismograms  to  constrain  the  velocity  and  density  structure  of  the  upper 
mantle.  Our  analysis  results  in  twelve  path-averaged  models  of  earth  structure  for  the  region 
of  the  Archaean  core.  Since  there  are  only  minor  differences  between  these  models  we  have 
averaged  them  to  form  a  1— D  earth  model  for  this  region.  We  then  compare  this  1— D  earth 
model  with:  (a)  velocity  and  density  estimates  derived  from  geochemical  studies  of  the  upper 
mantle  xenoliths  from  South  Africa,  (b)  previously  derived  upper  mantle  seismic  velocity  models 
for  southern  Africa,  and  (c)  seismic  velocity  models  derived  for  other  shield  regions.  The  goals  of 
our  study  are  to  place  better  constraints  on  the  velocity  structure  of  the  upper  mantle  beneath 
the  Archaean  core  of  southern  Africa;  to  determine  whether  the  predictions  of  in  situ  conditions 
based  on  the  geochemical  and  petrological  analysis  of  the  xenoliths  agree  with  seismic  measure¬ 
ments  in  the  same  region;  and  to  determine  if  the  xenolith  estimates  are  representative  of  the 
wider  region  of  the  Archaean  core.  In  the  main  text,  we  briefly  discuss  the  data  and  analytical 
method  we  follow,  summarize  the  seismological  results,  and  discuss  their  significance  in  light  of 
the  geochemical  analysis  of  nodules  from  the  kimberlites.  The  details  of  the  seismological  analysis 
and  the  geochemical  calculations  are  given  in  the  appendices. 

2.  GEOLOGICAL  AND  GEOPHYSICAL  CHARACTERISTICS  OF  SOUTH¬ 
ERN  AFRICA 

Southern  Africa  consists  of  a  well-defined  Archaean  core  surrounded  by  Proterozoic  mobile 
belts,  which  are  in  turn  surrounded  by  late  Proterozoic/early  Phanerozoic  mobile  belts  (Fig.  1). 
The  Archaean  core  consisting  of  the  Kaapvaal  and  Zimbabwe  cratons  and  the  Limpopo  mobile 
belt  has  formed  a  stable  unit  for  the  past  2.3  Ga  years  (McElhinny  &.  McWilliams,  1977).  There 
have  been  no  major  tectonic  events  within  southern  Africa  in  at  least  the  last  2  Ga  years  except 
in  the  Cape  Fold  belt  (Clifford,  1970).  The  geologic  and  tectonic  evolution  of  southern  Africa  is 
discussed  in  a  number  of  past  publications  (e.g..  Kroner,  1977;  Tankard  et  al,  1982;  de  Wit  et 
al,  1992). 

Two  observations  suggest  possible  differences  in  the  lithospheric  structure  within  southern 
Africa.  Although  kimberlites  are  found  over  much  of  this  region  (Pasteris,  1983),  diamond¬ 
bearing  kimberlites  have  a  more  restricted  distribution  and  occur  almost  totally  within  the  cratons 
(Gurney,  1990).  Since  the  diamond  stability  field  is  at  higher  pressure  (i.e.,  greater  depths) 
than  that  of  the  graphite  found  in  many  of  the  kimberlites,  the  distribution  of  diamond-bearing 
kimberlites  suggests  that  the  lithosphere  of  the  central  core  of  the  cratons  is  thicker  than  the 
lithosphere  of  the  surrounding  regions.  Some  South  African  kimberlites  are  as  old  as  1.75  Ga  years 
but  the  age  of  most  southern  African  kimberlites  clusters  around  80-90  million  years  (Dawson, 
1980).  However,  many  of  the  diamonds  are  Archaean  in  age  (Richardson  et  al,  1984),  implying 
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that  the  possible  differences  in  lithospheric  thickness  within  southern  Africa  suggested  by  the 
diamond-bearing  kimberlite  distribution  may  have  existed  since  the  early  stages  of  continental 
lithospheric  evolution.  This  in  turn  implies  that  these  deeper  portions  of  the  lithosphere  are  less 
dense  than  the  surrounding  mantle  since  they  have  not  delaminated  as  they  have  cooled. 

Heat  flow  observations  also  suggest  possible  differences  in  the  lithospheric  structure  within 
southern  Africa.  A  large  number  of  heat  flow  measurements  have  now  been  made  in  southern 
Africa.  Ballard  Pollack  (1987)  point  out  that  the  heat  flow  in  the  interior  of  the  Archaean  core 
is  typically  about  40  mWm”^  but  increases  to  about  60  mWm“^  at  the  craton  boundaries  with 
the  surrounding  Proterozoic  and  Pan-African  mobile  belts  and  as  high  as  70  mWm“^  within  the 
mobile  belts.  They  suggest  two  possible  mechanisms  to  explain  how  this  heat  flow  contrast  might 
arise:  a  shallow  geochemical  mechanism  resulting  from  a  difference  in  crustal  heat  production 
within  the  Archaean  core  and  surrounding  mobile  belts,  and  a  deeper  geodynamical  mechanism 
arising  from  a  thicker  lithospheric  root  beneath  the  Archaean  core  compared  to  that  beneath 
the  mobile  belts.  They  conclude  that  at  least  50%  and  possibly  100%  of  the  observed  heat  flow 
difference  results  from  a  thicker  lithospheric  root  beneath  the  Archaean  core  which  diverts  mantle 
heat  away  from  this  region  into  the  surrounding  mobile  belts.  To  achieve  this  they  suggest  that 
the  cratonic  root  beneath  southern  Africa  extends  to  a  depth  of  200  to  400  km.  Such  a  thick, 
low  density  lithospheric  root  for  the  craton  should  have  a  distinctive  seismic  signature. 

3.  DETERMINATION  OF  THE  VELOCITY  AND  DENSITY  STRUCTURE 
OF  THE  SOUTHERN  AFRICAN  CRUST  AND  UPPER  MANTLE 

3.1  Data  and  Method  of  Analysis 

We  have  used  seismograms  of  teleseismic  and  regional  earthquakes  (Table  1  &  2)  recorded 
at  six  digital  seismographs  (Table  3)  in  southern  Africa  to  determine  the  crust  and  upper  mantle 
structure  of  this  region.  The  locations  of  the  regional  earthquakes  and  seismograph  stations  are 
shown  in  Figure  2. 

We  employ  a  variety  of  seismological  techniques  to  determine  models  for  the  crust  and  upper 
mantle  velocity  structure.  We  determine  the  gross  features  of  the  crustal  structure  beneath  the 
seismograph  stations  using  teleseismic  receiver  function  modeling  (Owens  et  al.  1984).  The 
direct  teleseismic  P-wave  interacts  with  velocity  interfaces  and  gradients  in  the  crust  beneath  the 
seismograph  to  produce  a  seismogram  consisting  of  the  direct  P-wave  plus  P-wave  reverberations 
and  P-to-S  converted  phases  plus  their  reverberations.  Source  effects  can  be  largely  removed 
from  the  radial  component  by  deconvolving  from  it  the  vertical  component  seismogram  (Langston, 
1979).  For  laterally  homogeneous  structures  the  resulting  radial  component  seismogram,  called 
the  receiver  function,  primarily  contains  P-to-S  conversions  and  their  reverberations,  and  the 
resulting  tangential  component  is  zero.  Thus  the  ratio  of  the  tangential  to  radial  amplitude 
gives  an  indication  of  the  degree  of  scattered  energy  due  to  lateral  variation  in  earth  structure 
in  the  vicinity  of  the  seismograph.  We  have  modeled  the  receiver  functions  to  verify  that  crustal 
thickness  estimates  from  this  technique  do  not  differ  significantly  from  crustal  thickness  estimates 
from  nearby  seismic  refraction  data  where  they  exist  and  to  constrain  crustal  thickness  in  the 
vicinity  of  the  seismograph  stations  where  nearby  refraction  results  do  not  exist. 

We  determine  the  path-averaged  shear  wave  velocity  structure  for  the  crust  and  upper  man¬ 
tle  by  inverting  multi-mode  surface  waveforms  for  regional  earthquakes  using  the  technique  of 
Gomberg  Masters  (1988).  The  seismogram  is  represented  as  a  sum  of  locked-mode  travel¬ 
ing  waves  (surface  waves)  in  a  plane-layered  structure.  The  earth-flattening  transformations  of 
Biswas  Knopoff  (1970)  and  Biswas  (1972)  are  used  for  the  Love  and  Rayleigh  waves,  respec¬ 
tively,  and  physical  dispersion  is  accounted  for  using  the  correction  of  Liu  et  al.  (1976).  The 
difference  between  an  observed  and  synthetic  seismogram  can  be  directly  related  to  perturbations 
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in  the  model  velocity  and  density  structure  through  the  partial  derivatives  of  the  phase  velocity 
with  respect  to  the  model  parameters.  A  new  model  is  determined  and  the  process  repeated  until 
a  satisfactory  variance  reduction  is  achieved.  This  procedure  is  superior  to  inverting  surface  wave 
dispersion  curves  for  earth  structure  because  fundamental  mode  surface  waves,  whose  dispersion 
curves  can  readily  be  obtained,  on  their  own  are  not  very  sensitive  to  the  details  of  the  upper 
mantle  velocity  structure.  Dispersion  curves  of  higher  modes  are  more  sensitive  to  upper  mantle 
structure  but  are  difficult  to  determine  because  of  the  complex  interference  of  the  waves  arriving 
between  S  and  the  fundamental  mode  wavetrain.  By  directly  modeling  the  waveforms  this  inter¬ 
ference  is  taken  into  account  and  there  is  no  need  to  determine  spectral  properties  of  individual 
modes.  We  have  inverted  the  data  for  each  path  separately  but  we  have  simultaneously  inverted 
the  seismograms  from  the  vertical,  radial,  and  transverse  components. 

To  constrain  the  P-wave  velocity  gradient  in  the  upper  mantle  we  forward  model  Pn;  wave¬ 
forms  for  the  regional  earthquakes.  Pn/  is  a  combination  of  two  phases;  the  P^  phase  resulting 
from  P-wave  energy  turning  beneath  the  Moho  and  the  following  P;  mode  trapped  in  the  crustal 
wave  guide  (Helmberger  &  Engen,  1980).  P;  is  the  near  monochromatic  long-period  wavetrain 
often  seen  in  the  epicentral  distance  range  5-20°  following  the  Pn  phase  from  crustal  earthquakes 
(e.g.  Oliver  Si  Major,  1960;  Oliver,  1964).  P/  propagates  in  the  crustal  waveguide  as  a  partially 
trapped  P-Sv  wave  whose  phase  velocity  is  bracketed  by  the  compressional  and  shear  wave  ve¬ 
locity  of  the  upper  mantle  (Vs  <  cp,  <  Vp).  The  P-wave  is  therefore  post-critical  and  energy 
leaks  to  the  mantle  as  Sv,  resulting  in  a  decay  of  P/  with  range.  The  P/  waveform  is  primarily 
controlled  by  the  P-wave  travel  time  through  the  crust  (Shaw  Si  Orcutt,  1984)  resulting  in  a 
particular  interference  pattern  for  the  trapped  energy  and  a  characteristic  period  for  P;.  Because 
P/  is  only  weakly  affected  by  the  upper  mantle  velocity  structure  (Shaw  Si  Orcutt,  1984),  it 
forms  a  good  reference  for  comparing  the  amplitude  of  P„  or  direct  P  at  regional  distance  ranges 
(Clouser  Si  Langston,  1990),  which  is  very  sensitive  to  the  compressional  velocity  gradient  and 
attenuation  in  the  mantle  lid  (Hill,  1971;  1973).  We  forward  model  the  observed  Pn/  waveforms 
using  reflectivity  synthetics  (Fuchs  Si  Muller,  1971)  to  provide  a  constraint  on  the  upper  mantle 
compressional  velocity  gradients. 

3.2  Velocity  Structure  of  the  Southern  African  Crust 

This  study  is  primarily  concerned  with  the  average  velocity  structure  of  the  upper  mantle 
and  not  with  the  details  of  the  velocity  structure  of  the  crust.  However,  a  reasonably  accurate 
knowledge  of  crustal  structure  is  important  to  our  study  because  if  large  scale  lateral  variations 
in  crustal  structure  exist,  they  may,  if  not  properly  accounted  for,  map  into  the  structure  of  the 
upper  mantle.  A  number  of  seismic  refraction/reflection  studies  have  been  conducted  in  southern 
Africa  over  the  past  50  years.  Many  of  these  studies  have  used  mine  tremors  as  the  seismic  energy 
source.  These  have  the  advantage  that  they  are  a  rich  source  of  shear  wave  energy  but  have 
the  disadvantage  that  their  origin  times  and  locations  are  poorly  known  compared  to  those  of 
explosions.  Both  the  quality  of  these  data  and  the  methods  of  interpretation  have  varied  between 
these  studies,  so  it  is  not  easy  to  directly  compare  results  for  the  details  of  the  crustal  features. 
However,  the  crustal  thickness  determined  in  these  studies  can  be  compared.  All  of  the  cratonic 
seismic  profiles  are  for  the  Kaapvaal  Craton.  Durrheim  Si  Green  (1992)  analyzed  refraction 
data  for  two  seismic  profiles  in  the  central  Kaapvaal  Craton  by  matching  observed  and  synthetic 
seismograms,  and  their  study  provides  the  best  constraints  on  the  velocity  gradients  and  layering 
in  the  Kaapvaal  Craton  crust.  The  crust  in  this  region  has  the  following  features:  a  shallow  crust 
(<5  km  depth)  consisting  of  a  stratified  section  varying  in  compressional  velocity  from  3.5  to 
6.8  kms“^,  a  relatively  uniform  upper  crust  (5-12  km  depth)  with  velocities  ranging  from  6.0  to 

6.2  kms“^,  a  transition  from  upper  to  lower  crust  marked  by  a  change  in  velocity  gradient  at  12 
to  14  km  depth;  a  uniform  lower  crust  with  velocities  varying  from  about  6.4  kms"^  at  14  km 
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depth  to  about  6.7  kms~^  at  about  32  km  depth,  and  a  gradational  crust/mantle  transition  near 
35  km  depth  with  velocities  increasing  from  less  than  7  kms“^  to  greater  than  8  kms“^  over  a 
2  to  3  km  depth  range.  The  crustal  shear  wave  velocity  from  traveltime  studies  varies  from  3.66 
kms-i  (Gane  et  al,  1956)  to  4.15  kms-^  (Hales  &  Sacks,  1959).  Durrheim  &  Green  (1992) 
find  the  average  shear  wave  velocity  of  the  crust  to  be  3.8  kms-T  The  Pn  velocity  beneath  the 
Kaapvaal  Craton  is  8.0  -  8.2  kins-^  (Willmore  et  al,  1952;  Hales  &  Sacks,  1959;  Durrheim  & 
Green,  1992).  Willmore  et  al  (1952)  observe  to  have  a  velocity  of  4.83  kms-\  significantly 
higher  than  the  4.70-4.73  km s“^observed  by  Gane  et  al  (1956)  and  Durrheim  &.  Green  (1992). 

To  determine  the  crustal  structure  in  the  immediatevicinity  of  thefour  broadband  seismograph 
stations,  we  have  analyzed  teleseismic  receiver  functions.  This  is  particularly  important  in  the 
case  of  LSZ  (Fig.  2)  where  no  nearby  refraction  data  exist  and  the  crustal  structure  is  otherwise 
unknown.  Appendix  A  summarizes  the  receiver  function  analysis.  In  analyzing  the  receiver 
function  data  we  are  primarily  concerned  with  determining  the  crustal  thickness  and  average 
crustal  velocity.  Since  the  receiver  functions  consist  of  few  events  we  have  concentrated  on 
fitting  the  early  arrivals  and  have  not  emphasized  fitting  later  weak  arrivals. 

Figure  3  shows  the  stacked  radial  and  transverse  receiver  functions  for  BOSA,  along  with  the 
i  1  <T  bounds  determined  from  the  variance  of  the  stack.  These  bounds  provide  a  useful  measure 
of  the  coherence  of  the  arrivals  comprising  the  stack.  The  radial  receiver  function  shows  a  large 
amplitude  arrival  4  seconds  after  the  main  direct  P  arrival.  There  are  two  smaller  amplitude 
arrivals;  a  positive  arrival  ~  14  seconds  after  direct  P  and  a  negative  arrival  ~  18  seconds  after 
direct  P.  The  amplitudes  of  these  arrivals  on  the  radial  are  large  compared  with  the  amplitudes 
of  the  tangential  at  the  same  time,  if  we  interpret  the  first  of  these  later  arrivals  as  a  P-to-S 
conversion  at  the  Moho  (Ps)  ,  assume  an  average  crustal  P-wave  velocity  of  6.4  kms~^  and  a 
Poisson’s  ratio  of  0.25,  the  timing  of  this  arrival  suggests  a  35  km  thick  crust  beneath  BOSA. 

We  have  inverted  the  BOSA  receiver  function  data  in  Figure  3  and  similar  data  for  LBTB, 
LSZ,  and  SUR  using  the  linearized  inversion  procedure  of  Ammon  et  al  (1990).  The  receiver 
function  result  for  BOSA  is  shown  in  the  inset  in  Figure  3  and  the  results  for  the  other  three  sites 
are  shown  in  Figures  A2-5.  The  receiver  function  crustal  model  for  BOSA  has  a  steep  gradient 
in  the  shallow  crust  with  the  shear  wave  velocity  increasing  from  2.35  kms~^  at  the  surface  to 
3.67  kms  ^  at  2.6  km  depth.  The  crystalline  crust  is  composed  of  a  26  km  thick  section  with 
an  average  velocity  of  3.67  kms"^  There  is  a  steep  gradient  between  29  and  35  km  depth  with 
the  Moho  transition  between  velocities  of  3.98  and  4.75  kms~^  occurring  from  35  to  38  km 
depth.  The  synthetic  receiver  function  fits  the  general  features  of  the  observation  but  not  the 
details.  Much  better  fits  were  achieved  than  that  shown  in  Figure  3  by  allowing  more  complex 
earth  models.  However,  since  the  receiver  functions  were  composed  of  few  events  and  because 
our  purpose  was  to  confirm  the  gross  crustal  features  and  not  to  examine  the  details  of  the  crustal 
velocity  structure,  we  have  chosen  only  to  model  the  main  features  of  the  receiver  functions. 

Four  seismic  refraction  studies  (Willmore  et  al,  1952;  Gane  et  al,  1956;  Hales  ^  Sacks,  1959; 
Durrheim  &  Green,  1992)  show  the  crust  is  36il  km  thick  over  a  wide  region  of  the  central 
part  of  the  Kaapvaal  Craton  (Fig.  4).  The  teleseismic  receiver  function  result  for  BOSA  in  the 
SW  part  of  the  craton  shows  the  crust  is  35  km  thick  in  this  region  also.  The  receiver  function 
result  for  LBTB  (Fig.  A2)  in  the  NW  part  of  the  craton  shows  the  crust  in  this  region  is  44  km 
thick,  significantly  thicker  than  to  the  south.  This  difference  is  clearly  seen  from  comparing  the 
Ps-P  time  differences  for  the  BOSA  and  LBTB  receiver  functions;  this  time  difference  is  nearly 
one  second  greater  for  LBTB  (Fig.  A2)  than  for  BOSA  (Fig.  3).  Wright  &  Hall  (1990)  suggest 
that  their  seismic  reflection  data  from  the  Kalahari  of  western  Botswana  to  the  NW  of  LBTB 
show  the  crust  thickening  to  the  SE  approaching  the  Kaapvaal  Craton.  They  observe  a  14  sec 
two-way  travel  time  from  a  reflector  they  identify  as  the  Moho.  If  the  average  crustal  P-wave 
velocity  is  6.4  kms  ^  in  this  region  as  it  is  to  the  south,  this  corresponds  to  a  Moho  at  about 
45  km  depth.  Therefore  it  may  be  that  beneath  much  of  the  Kaapvaal  Craton  the  crust  is  36±1 
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km  thick  but  is  significantly  thicker  in  the  I\IW  part  of  the  craton. 

The  crustal  structure  of  the  Zimbabwe  Craton  and  the  Limpopo  mobile  belt  is  poorly  known. 
The  little  seismic  refraction  data  that  exist  suggest  that  the  crust  beneath  the  Zimbabwe  Craton 
is  about  40  km  thick  and  that  the  crust  beneath  the  Limpopo  mobile  belt  is  34-36  km  thick 
(Stuart  &  Zengeni,  1987;  Durrheim  et  ah,  1992).  These  values  agree  with  the  crustal  thickness 
estimates  from  gravity  and  electrical  studies  (de  Beer  &  Stettler,  1992:  Gwavava  et  al,  1992). 
Beneath  the  crust  of  the  Limpopo  there  exists  a  5-10  km  thick,  8.1  kms"^  layer  overlying  a  14 
km  thick  layer  where  Vp  drops  to  7.6  kms"^  (Durrheim  et  al,  1992);  beneath  this  there  is  a 
positive  velocity  gradient  with  velocities  reaching  8.5  km  8“^  at  a  depth  of  110  km. 

Seismic  refraction/reflection  studies  show  that  the  crust  beneath  the  Proterozoic/early 
Phanerozoic  mobile  belts  is  thicker  (42-47  km)  than  beneath  the  the  Archaean  core  (Baier 
et  al,  1983;  Green  &.  Durrheim,  1990;  Wright  &  Hall,  1990).  In  addition,  there  is  a  suggestion 
(Durrheim  &.  Mooney,  1991)  that  the  lower  crust  beneath  the  mobile  belts  contains  a  thick,  high 
velocity  (Vp  >  7.0  kms"^)  section  not  observed  in  the  lower  crust  of  the  cratons.  However, 
most  of  the  analysis  of  seismic  refraction  data  in  southern  Africa  predate  the  use  of  synthetic 
seismogram  modeling  and  hence  this  proposed  difference  in  the  lower  crust  between  the  mobile 
belts  and  the  cratons  requires  substantiation.  SUR  is  near  the  boundary  of  the  Namaqua  Province 
and  the  Cape  Fold  Belt.  Inversion  of  two  receiver  function  data  (Fig.  A3-4)  shows  that  the  crust 
beneath  this  site  is  52il  km  thick.  Green  Durrheim  (1990)  find  the  crustal  thickness  from 
seismic  refraction  measurements  in  the  Namaqua  metamorphic  complex  300  km  NW  of  SUR 
is  significantly  thinner  (42  km).  LSZ  is  located  north  of  the  Zimbabwe  Craton  in  the  Irumide 
mobile  belt.  We  have  inverted  the  single  event  receiver  function  from  LSZ  (Fig.  A5)  and  the 
result  suggests  a  crustal  thickness  of  about  42  km.  Since  the  receiver  function  consists  of  only 
a  single  event,  we  have  less  confidence  in  this  result  than  in  the  other  receiver  function  results; 
however,  this  is  the  only  crustal  thickness  estimate  available  for  this  region. 

Figure  4  summarizes  the  data  for  crustal  thickness  in  southern  Africa.  From  these  data  we 
form  an  average  1-D  velocity  model  for  the  southern  African  crust  to  use  as  a  starting  model 
for  the  regional  waveform  inversion.  This  crustal  model  is  41.4  km  thick  and  has  the  following 
main  features:  a  1.4  km  thick  sedimentary  layer  with  P-wave  velocity  4.2  kms”^,  a  5  km  thick 
upper  crust  with  P-wave  velocity  6.2  kms-\  a  20  km  thick  mid-crust  with  P-wave  velocity 
6.4  kms■^  and  a  15  km  thick  lower  crust  with  P-wave  velocity  6.8  kms“T  Poisson’s  ratio  is 
assumed  to  be  0.25  throughout  the  crust. 

3.3  Velocity  Structure  of  the  Southern  African  Mantle 

3.3.1  Body-wave  constraints 

Knowledge  of  the  compressional  and  shear  wave  speeds  directly  beneath  the  Moho  (Vp„  and 
Vs„)  provide  important  constraints  in  the  inversion  of  the  regional  surface  waveforms.  More 
precise  knowledge  of  Vs„  beneath  the  Archaean  core  is  particularly  important  in  our  study  since 
surface  wave  propagation  is  primarily  controlled  by  the  shear  wave  structure.  Both  the  fact  that 
Sn  travels  with  a  constant  surface  velocity  near  that  of  the  upper  mantle  shear  speed  and  that 
Sn  has  a  substantial  high  frequency  component  suggest  that  it  propagates  as  S-wave  energy 
trapped  in  the  upper  mantle  lid,  and  is  thus  diagnostic  of  the  upper  mantle  shear  wave  speed. 
Independent  knowledge  of  the  upper  mantle  shear  wave  speed  is  important  since  it  can  aid  in 
distinguishing  between  velocity  models  with  low  average  upper  mantle  shear  wave  velocities  but 
no  LVZ  and  velocity  models  containing  a  higher  velocity  upper  mantle  lid  above  a  shear  wave 
LVZ,  both  of  which  could  agree  with  the  surface  wave  observations. 

Vp„  measurements  from  seismic  refraction  profiles  in  several  areas  of  southern  Africa  (Fig.  4) 
show  Vp„  is  8.1  ±  0.1  kms-T  The  measured  S„  velocity  observed  along  three  of  the  refraction 
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profiles  in  the  Kaapvaal  Craton  (Willmore  et  al,  1952;  Gane  et  al,  1956;  Durrheim  ^  Green, 
1992)  varies  from  4.70  to  4.83  kms  Gumper  ^  Pomeroy  (1970)  measured  Vp^,  8.06  kms“^ 
and  Vs„  4.55-4.72  knas"^  using  earthquakes  from  over  a  wider  region  of  southern  and  eastern 
Africa.  However,  several  of  the  paths  they  examined  traversed  the  East  African  Rift.  Molnar  L 
Oliver  (1969)  and  Gumper  &  Pomeroy  (1970)  note  that  is  only  observed  for  paths  in  southern 
Africa  and  for  paths  south  of  10°S  crossing  the  rift.  The  observed  values  of  Vs„  (Gumper  & 
Pomeroy,  1970)  over  this  wider  region  of  southern  Africa  span  the  range  observed  for  shields 
(e.g.  Brune  &:  Dorman,  1963;  The  Canadian  Shield)  to  regions  of  active  continental  rifting  (e.g. 
Priestley  &  Brune  1978;  The  Basin  and  Range). 

To  improve  the  constraint  on  Vs„  in  southern  Africa  we  constructed  a  composite  record  section 
of  seismograms  from  southern  African  events  recorded  in  the  distance  range  300-2,000  km.  The 
propagation  paths  for  these  events  are  primarily  confined  to  the  stable  regions  of  southern  Africa 
(Fig.  2  &  Table  1)  although  some  of  the  events  occur  in  the  southernmost  part  of  the  rift. 
These  seismograms  (Fig.  5)  show  a  clear  S„  phase  with  a  travel  time  curve  given  by  T=12.0 
+  ^/4-64.  Vs„  ~  4.64  km  s“^  falls  in  the  middle  of  the  Gumper  &  Pomeroy  (1970)  Vs„ 
range  But  is  low  compared  with  Vs„  observed  for  shorter  range  refraction  recordings  within  the 
Kaapvaal  Craton  (Willmore  et  al,  1952;  Gane  et  al.  1956;  Durrheim  &  Green,  1992)  and  for 
most  shields  (Huestis  et  al,  1973).  In  the  surface  wave  inversion  we  use  Vs  =  4.64  km 
as  a  constraint  on  the  S-wave  velocity  immediately  below  the  Moho.  In  modeling  events  with 
most  of  their  path  length  within  the  Archaean  core  or  mobile  belts  of  southern  Africa  we  require 
Vs„  >  4.64  kms“^  but  relaxed  this  requirement  for  events  with  larger  fractions  of  their  path 
length  in  the  southern  portions  of  the  East  African  Rift. 

3.3.2  Modeling  of  Regional  Surface  Waves 

We  determine  the  mantle  shear  velocity  structure  using  the  differential  seismogram  technique 
of  Gomberg  &.  Masters  (1988).  We  have  inverted  twelve,  three-component  long  period  seismo¬ 
grams  from  eight  regional  earthquakes  (<  30°)  in  southern  Africa  (Fig.  2  &  Table  1)  recorded 
at  the  digital  stations  BOSA,  LBTB,  LSZ.  SLR,  and  SUR  (Table  3).  This  has  resulted  in  twelve 
path-averaged  earth  models.  The  propagation  paths  are  primarily  confined  to  the  cratons  and 
portions  of  the  surrounding  mobile  belts;  however,  some  paths  also  sample  small  portions  of  the 
east  African  rift.  Appendix  B  describes  in  more  detail  the  analysis  and  summarizes  the  inversion 
of  the  regional  surface  waveforms. 

The  seismogram  of  the  July  18,  1986  (Event  2,  Table  1  &  2)  earthquake  recorded  at  SLR  is 
typical  of  the  data  we  have  inverted  for  earth  structure  and  we  use  it  as  an  example  of  our  analysis 
procedure.  This  event  occurred  on  the  northern  boundary  of  the  Zimbabwe  Craton;  thus,  the 
propagation  path  to  SLR  was  confined  almost  entirely  to  the  Archaean  core.  The  event  is  almost 
due  north  of  SLR  and  the  seismograms  (Fig.  6)  are  therefore  close  to  being  naturally  rotated 
(N-S  radial;  E-W  tangential).  The  short-period  P-wave  consists  of  an  impulsive  (group 
velocity  7.80  kms-^)  and  no  significant  Pg  phase.  This  observation  is  typical  of  short-period 
P-wave  propagation  in  stable  regions  (Evernden,  1967;  Langston,  1982).  The  S„  phase  (group 
velocity  4.46  kms  ^)  is  followed  by  a  large  amplitude  Lg  phase  (beginning  group  velocity  ~  3.67 
kms  ).  The  long-period  seismograms  are  relatively  simple.  The  first  arrival  on  the  vertical 
and  N-S  (radial)  primarily  consists  of  the  wave.  The  surface  waves  begin  with  the  higher 
Rayleigh  and  Love  modes  near  the  predicted  S-wave  arrival  time  (~  235  sec)  and  culminate  with 
the  fundamental  Love  mode  Airy  phase  at  about  300  seconds  and  fundamental  Rayleigh  mode 
Airy  phase  at  about  335  seconds. 

The  starting  earth  model  for  the  waveform  inversion  consists  of  the  average  crustal  model 
discussed  in  section  3.2  and  the  upper  mantle  S-wave  model  from  Bloch  et  al  (1969)  (Fig. 
7c).  The  upper  mantle  P-wave  velocity  was  related  to  the  S-wave  velocity  assuming  a  Poisson's 
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ratio  of  0.25  and  the  upper  mantle  density  was  related  to  the  P-wave  velocity  assuming  p  = 
0.77  +  0.32V^.  Figure  7a  compares  the  long-period  seismograms  with  synthetic  seismograms 
(fundamental  and  eight  higher  modes)  computed  for  Event  2  at  SLR  using  the  ISC  location 
(Table  1)  and  the  CMT  source  parameters  (Table  2)  and  the  starting  earth  model.  The  phase  of 
the  synthetic  higher  mode  matches  the  observed  higher  mode  although  the  synthetic  amplitude 
is  larger  than  the  observed  amplitude.  The  long  period  synthetic  fundamental  mode  phase  and 
amplitude  are  in  reasonable  agreement  with  the  observed  waveform  but  the  fundamental  mode 
Airy  phase  is  in  poor  agreement.  Before  inverting  the  surface  waveforms  for  a  new  earth  model  we 
first  checked  to  verify  that  the  observed  misfit  could  not  simply  be  explained  by  reasonable  errors 
in  the  hypocentral  location  or  mechanism  for  the  event.  The  ISC  location  accuracy  in  southern 
Africa  is  uncertain  but  a  comparison  of  the  teleseismic  locations  of  the  ISC  with  regional  locations 
of  the  South  African  Geological  Survey  for  large  rock  bursts  suggest  that  the  ISC  location  errors 
do  not  greatly  exceed  ±  10  km  for  the  three  hypocentral  coordinates.  Altering  the  hypocenter  by 
these  amounts  and  the  source  parameters  (strike,  slip,  dip)  ±10°  from  the  CMT  values  did  not 
produce  a  better  fit  between  the  observed  and  synthetic  seismograms  computed  for  the  starting 
model. 

The  three-component  long-period  seismograms  were  inverted  for  the  path-average  earth 
model  using  the  waveform  fitting  approach  of  Gomberg  &.  Masters  (1988).  The  new  earth  model 
is  shown  in  Figure  7c.  The  mantle  consists  of  a  76  km  thick  upper  mantle  lid  with  Vs„  = 
4.64  kms“^  directly  beneath  the  Moho.  The  two  layer  lid  corresponds  to  a  shear  wave  velocity 
gradient  of  about  0.0011  s"^.  Below  the  lid  is  a  substantial  upper  mantle  LVZ.  At  the  bottom 
of  the  lid  Vs  equals  4.72  kms~\  but  drops  to  4.10  kms"^  at  200-330  km  depth.  The  synthetic 
seismograms  computed  for  the  inversion  earth  model  are  shown  in  Figure  7b.  We  found  that  the 
fit  was  improved  by  reducing  the  hypocentral  depth  of  the  earthquake  from  17  to  12  km.  The 
phase  of  the  higher  mode  wave  train  arriving  between  about  240  and  300  sec  is  well  fit  by  the 
synthetic  seismogram  although  the  synthetic  amplitude  is  still  somewhat  high.  The  low  frequency 
portion  of  the  fundamental  mode  wavetrain  is  well  fit  and  the  phase  of  the  fundamental  mode 
Airy  phase  is  matched;  however,  the  amplitude  of  the  main  peak  on  the  vertical  component  is 
underestimated  by  about  40%. 

In  Appendix  B  we  discuss  a  series  of  tests  we  have  conducted  to  assess  the  effects  of  errors  in 
our  assumptions  and  in  the  data  on  the  inversion  earth  model.  We  assume  that  the  hypocenter 
may  be  in  error  by  ±  10  km  in  distance  and  depth  and  that  the  source  mechanism  may  be  in 
error  by  ±10°  in  strike,  dip,  and  slip.  We  then  reinvert  the  surface  waveforms  to  find  the  effects 
of  these  values  on  the  earth  structure.  The  results  of  these  tests  are  shown  as  the  bounds  on  the 
velocities  and  density  in  Figure  7c.  We  have  also  conducted  tests  to  evaluate  the  resolution  of 
the  data  in  controlling  the  thickness  of  the  upper  mantle  lid  and  the  existence  of  a  LVZ.  While 
the  details  of  the  LVZ  velocity  structure  varied,  we  were  unable  to  find  shear  wave  upper  mantle 
models  which  satisfied  the  Vs„  data  from  Figure  5  and  which  did  not  contain  a  S-wave  LVZ.  We 
therefore  conclude  that  at  least  an  upper  mantle  shear  wave  LVZ  exists  beneath  the  Archaean 
core  of  southern  Africa. 

3.3.3  Modeling  of  Pn/  Waveforms 

Surface  wave  observations  do  not  constrain  upper  mantle  compressional  wave  velocity  struc¬ 
ture.  We  determine  the  compressional  velocity  structure  of  the  upper  mantle  by  forward  modeling 
Pn/  waveforms  from  the  seismograms  of  the  same  eight  regional  earthquakes  used  in  the  surface 
waveform  inversion.  We  follow  an  approach  similar  to  Clouser  <Si  Langston  (1990).  Appendix  B 
also  summarizes  our  forward  modeling  of  the  Pn/  waveform  data. 

Figure  8  is  an  enlargement  of  the  P-wave  for  Event  2  recorded  at  SLR.  To  gain  insight 
into  the  composition  of  the  P-wave  we  first  smoothed  the  coarse  “blocky”  inversion  model  by 


increasing  the  number  of  layers  and  changing  some  of  the  velocity  steps  in  the  coarse  model  into 
gradients.  We  then  computed  generalized  ray  synthetic  seismograms  for  this  smoothed  surface 
wave  inversion  model  for  the  separate  arrivals  making  up  the  Event  2  P-wave.  Based  on  these 
synthetics  the  first  arrival  is  the  direct  P-wave  and  the  second  arrival  is  composed  of  upper 
mantle  turning  rays  from  the  S-to-P  conversion  at  the  Moho  beneath  the  source  and  from  the 
depth  phases  pP  and  sP.  Comparing  the  individual  amplitudes  of  these  phases  shows  that  sP 
is  primarily  responsible  for  the  observed  waveform  for  this  event.  Figure  9a  shows  reflectivity 
synthetics  computed  for  the  smoothed  surface  wave  inversion  model  compared  with  the  Event  2 
long-period  SLR  seismograms.  Whereas  the  overall  fit  of  the  Pn/  waveform  is  good,  the  amplitude 
of  the  synthetic  first  arrival  is  small  compared  to  that  of  the  observed  first  arrival  suggesting  that 
there  is  insufficient  energy  present  in  the  synthetic  mantle  turning  wave  at  this  range. 

We  replaced  the  mantle  lid  P-wave  velocity  structure  derived  from  the  surface  wave  modeling 
with  the  mantle  lid  P-wave  structure  SACM06c  found  by  Clouser  &  Langston  (1990)  and  recom¬ 
puted  the  reflectivity  synthetic  seismograms.  This  produces  a  much  better  match  to  the  observed 
first  arrival.  We  then  used  a  trial-and-error  forward  modeling  approach  of  matching  the  synthetic 
waveform  to  the  observed  waveforms  by  adjusting  the  upper  mantle  P-wave  velocity  gradient. 
The  Event  2-SLR  path  is  1038  km  long,  similar  to  the  most  distant  event  studied  by  Clouser 

Langston  (1990).  The  match  between  the  synthetic  and  observed  seismograms  for  this  event 
imply  that  there  is  a  positive  compressional  velocity  gradient  in  the  upper  mantle  lid  beneath  the 
Archaean  core,  and  that  the  base  of  the  "P-wave  lid”  is  at  least  50  km  deeper  than  the  base 
of  the  “S-wave  lid”.  However  in  modeling  events  at  greater  distance  ranges  (see  Appendix  B) 
we  found  that  the  top  of  the  "P-wave”  LVZ  must  be  even  deeper  and  that  a  good  fit  to  all 
Pn/  data  was  obtained  for  the  case  of  no  P-wave  LVZ.  The  final  earth  model  for  the  analysis  of 
seismograms  from  event  2  is  depicted  in  Figure  9c  and  tabulated  in  Table  4.  The  smoothing  of 
the  upper  mantle  velocity  and  density  gradients  has  also  improved  the  amplitude  fit  of  the  higher 
mode  surface  waves  (Fig.  9b). 

Our  study  of  Pn/  waveforms  confirms  Clouser  Langston’s  (1990)  observation  that  a  positive 
P-wave  velocity  gradient  exists  in  the  upper  mantle  beneath  the  Archaean  core  but  we  find  the 
gradient  is  0.0008  s“^,  smaller  than  in  model  SACMOGc.  However,  our  data  imply  that  if  a 
P-wave  LVZ  exists,  it  is  much  deeper  than  Clouser  Si  Langston’s  (1990)  upper  bound  of  170  km. 

3.3.4  One-Dimensional  Earth  Model  for  the  Southern  African  Craton 

Data  for  eleven  other  source-receiver  paths  (Fig.  2)  have  been  analyzed  in  the  same  manner 
as  described  above,  and  the  results  for  these  are  given  in  Appendix  B.  The  synthetic  fits  to 
the  observed  waveforms  shown  above  for  the  July  18,  1986  event  are  typical;  some  synthetic 
fits  shown  in  Appendix  B  are  better  than  those  for  Event  2  while  others  are  somewhat  poorer. 
The  propagation  paths  for  these  events  primarily  sample  the  stable  region  of  southern  Africa  but 
waveforms  for  some  paths  sample  larger  portions  of  the  Proterozoic  belts  than  other  paths.  There 
are  differences  between  the  models  for  the  twelve  paths;  however,  these  are  in  details  and  not  in 
the  main  upper  mantle  features,  suggesting  that  there  are  no  major  lateral  variations  in  the  upper 
mantle  structure.  We  have  averaged  these  models  to  obtain  a  1-D  earth  model  for  the  Archaean 
core  of  southern  Africa  (Fig.  10  Si  Table  5).  This  earth  model  has  the  following  features.  The 
42  km  thick  crust  consists  of  four  layers;  a  2  km  thick  5.05  kms“^  surface  layer,  a  5  km  thick 
6.08  kms“^  upper  crust,  a  20  km  thick  6.30  kms“^  mid-crust,  and  a  15  km  thick  6.73  kms”^ 
lower  crust.  Below  the  crust  there  is  an  80  km  thick  upper  mantle  lid.  Vp„  and  Vs„  are  8.09 
and  4.62  kms“\  respectively,  and  the  compressional  and  shear  velocity  gradients  in  the  lid  are 
0.0008  s“^  and  0.0013  s"^,  respectively.  Below  the  lid  there  is  a  substantial  shear  wave  LVZ  with 
the  S-wave  velocity  dropping  to  4.32  km  s~^  at  250  km  depth.  Below  125  km  depth  the  P-wave 
gradient  increases  to  0.0015  and  increases  again  to  0.0035  s“^  between  250  km  depth  and 


the  410  discontinuity. 


4.  DISCUSSION 

4.1  Comparison  Seismically  Determined  Velocities  with  those  Estimated  from 
Kimberlite  Nodule  Data 


An  important  reason  for  investigating  seismic  wave  propagation  across  the  Archaean  cratons 
of  southern  Africa  is  that  a  great  deal  of  information  is  now  available  about  the  composition, 
pressure  and  temperature  at  depth  beneath  these  cratons  from  studies  of  mantle  nodules,  and  the 
kimberlitic  host  magmas  that  transported  them  to  the  surface.  The  pressure  and  temperature 
at  which  the  nodules  last  equilibrated  can  be  estimated  from  two  metamorphic  reactions.  The 
intersolubility  of  ortho-  and  clinopyroxene  is  principally  dependent  on  temperature,  whereas  the 
solubility  of  AI2O3  in  orthopyroxene,  in  the  presence  of  garnet,  is  principally  controlled  by  pressure. 
Finnerty  &  Boyd  (1987)  review  these  and  other  reactions  that  can  be  used  as  geothermometers 
and  geobarometers.  These  estimates  have  been  combined  with  thermal  models  of  the  lithosphere 
(McKenzie,  1989)  to  obtain  the  thicknesses  of  the  mechanical  and  thermal  boundary  layers 
beneath  the  Kaapvaal  Craton  and  surrounding  regions.  This  model  is  used  here,  with  a  mechanical 
boundary  layer  thickness  of  165  km  and  a  viscosity  of  4  x  10^®  in  the  thermal  boundary 

layer,  to  estimate  the  temperature  at  depth.  It  gives  a  thermal  boundary  layer  thickness  of  36 
km,  a  lithospheric  thickness  of  185  km  and  a  mantle  heat  flow  of  23.1  mWm“^. 

The  seismic  velocities  and  density  depend  on  the  composition  and  mineralogy  of  the  mantle, 
as  well  as  on  the  pressure  and  temperature.  The  nodules  provide  samples  of  mantle  minerals, 
and  the  average  composition  of  the  more  important  minerals  is  listed  in  Table  Cl  (see  McKenzie 
&  0  Nions,  1991,  for  details).  The  proportions  of  these  minerals  present  at  depth,  known  as  the 
modal  mineralogy,  can  be  determined  either  by  point  counting,  or  by  finding  the  proportions  of 
the  minerals  that  best  fit  the  bulk  composition  of  the  nodules.  The  bulk  composition  can  either  be 
measured  or  calculated  from  a  compositional  model  of  the  source.  Tainton  &.  McKenzie  (1994) 
showed  that  all  these  approaches  gave  similar  estimates  of  modal  mineralogy.  Unfortunately  the 
modal  mineralogy,  mineral  and  bulk  composition  have  not  all  been  reported  for  the  same  nodules. 
We  therefore  used  the  melting  models  of  Tainton  &.  McKenzie  (1994)  that  are  consistent  with 
all  the  information  available  from  the  nodules  and  their  host  kimberlites  to  estimate  the  modal 
mineralogy  as  a  function  of  depth.  We  used  the  models  for  four  nodule  suites:  PKP  (phogopite 
K— richterite  peridotite),  and  GP  and  GPP  (garnet  peridotite  and  garnet  phlogopite  peridotite) 
from  Erlank  et  al.  (1987),  and  the  depleted  and  fertile  nodules  from  Nixon  et  al.  (1981),  and 
obtained  the  modal  mineralogy  by  minimising  the  misfit  between  the  bulk  composition  and  that 
calculated  from  the  phases  in  Table  Cl.  The  minimization  used  6  phases,  olivine,  orthopyroxene, 
clinopyroxene,  garnet,  phlogopite  and  chrome-spinel,  and  the  misfit  function  was  the  sum  of  the 
squares  of  the  differences  between  the  oxide  compositions,  with  the  difference  for  K2O  weighted 
by  a  factor  of  ten.  The  modal  mineralogy  of  the  four  nodule  suites  is  listed  in  Table  Cl,  together 
with  the  parameters  used  to  calculate  the  density  and  the  seismic  velocities  from  the  modal 
mineralogy.  In  all  cases  the  ratio  of  MgO  to  FeO  was  fixed  to  be  the  same  as  that  in  the  nodule 
suite  in  olivine,  orthopyroxene  and  clinopyroxene.  The  method  used  to  obtain  p,  Vp  and  Vs  is 
outlined  in  the  Appendix  C,  and  can  be  used  to  calculate  these  parameters  for  any  rock  that 
consists  of  the  minerals  in  Table  Cl. 

The  nodule  compositions  from  southern  Africa  show  that  the  composition  of  the  upper  mantle 
sampled  by  the  kimberlites  varies  with  depth.  Unfortunately,  accurate  depth  estimates  cannot 
be  made  for  the  four  suites  listed  in  Table  Cl  because  their  mineral  compositions  have  not  been 
reported.  But  Waters  <&:  Erlank  (1988)  suggest  that  the  PKP  nodules  come  from  depths  of  about 
100  -  130  km.  Depleted,  GP  and  GPP  nodules  give  pressure  estimates  from  100  -  150  km  and 


fertile  nodules  from  160  -  200  km  (see  Tainton  &.  McKenzie,  1994). 

The  velocities  and  density  computed  from  the  nodules  are  given  in  Table  6;  the  details  of 
the  calculations  are  given  in  Appendix  C.  The  petrological  estimates  are  compared  with  the 
seismic  velocities  and  density  in  Figure  10  and  the  lid  values  are  in  excellent  agreement,  which 
gives  strong  support  for  the  proposed  upper  mantle  composition  and  also  implies  that  the  upper 
mantle  conditions  deduced  from  the  petrological  data  are  representative  of  the  Archaean  core  of 
southern  Africa, 

However,  we  predict  significantly  higher  velocities  and  densities  from  the  deepest  nodules 
(Table  6)  compared  to  the  velocities  and  density  at  corresponding  depths  in  our  seismic  model. 
The  largest  difference  is  for  the  S-wave  velocity.  The  depth  to  the  top  of  the  S-wave  LVZ  in 
our  model  for  southern  Africa  is  122  km  and  we  argue  in  Appendix  B  that  this  is  resolvable  to 
about  ±30  km.  Thus  even  the  maximum  depth  to  the  top  of  the  LVZ  is  in  disagreement  with 
the  deepest  nodule  data.  This  discrepancy  is  greatest  for  the  S-wave  velocity  and  the  reason 
for  this  may  be  due  to  the  presence  of  small  melt  fractions  just  below  the  seismically  defined 
upper-mantle  lid.  McKenzie  (1989)  has  shown  that  small  melt  fractions  are  likely  to  be  present 
beneath  the  lithosphere  and  these  will  separate  from  the  matrix  when  there  is  more  than  about 
10“^%  melt  present  if  their  viscosity  is  0.1  Pa  s.  Once  separated  these  melts  will  migrate  upwards 
but,  since  such  melts  movement  cannot  transport  heat,  they  will  solidify  within  the  mechanical 
boundary  layer  at  temperatures  between  950°  and  750°  (McKenzie,  1989).  The  corresponding 
depths  can  be  estimated  from  the  nodule  data  and  for  the  Kaapvaal  Craton  these  temperatures 
correspond  to  depths  of  120  to  100  km  if  the  dihedral  angle  is  as  small  as  zero.  Below  this  depth 
this  small  melt  fraction  can  coat  the  grain  boundaries  of  the  matrix  with  a  thin  film.  Such  a 
film  has  been  assumed  to  be  absent  when  the  bulk  elastic  parameters  are  calculated  from  the 
modal  mineralogy.  Though  its  presence  would  have  little  effect  on  the  density  and  bulk  modulus, 
the  same  is  not  true  of  the  shear  modulus.  Lubrication  of  the  grain  boundaries  could  have  a 
significant  effect  on  the  shear  modulus,  which  controls  Vs  (O’Connell  &  Budiansky,  1974). 

4.2  Geological  Implications  of  the  Velocity  Structure 

Though  southern  Africa  has  not  been  deformed  since  the  Archaean,  it  has  undergone  regional 
uplift  of  at  least  1  km  over  very  large  regions.  McKenzie  (1984)  and  Cox  (1993)  have  suggested 
that  this  uplift  results  from  the  underplating  of  Karoo  basalts.  This  suggestion  is  consistent  with 
our  velocity  model  if  such  material  is  emplaced  in  the  lower  crust.  The  velocities  and  density  of 
the  15  km  thick  layer,  at  depths  of  27-42  km,  are  typical  of  those  of  basalt,  and  emplacement  of 
such  a  layer  would  produce  an  uplift  of  1.4  km. 

The  velocity  of  the  S  waves  in  the  high  velocity  mantle  lid  agree  better  with  those  calculated 
for  the  PKP  nodules.  Though  the  difference  between  the  velocities  of  different  nodule  types  is 
scarcely  larger  than  the  likely  modelling  errors,  only  the  PKP  nodules  come  from  depths  between 
52  and  122  km,  corresponding  to  the  high  velocity  lid.  The  velocities  of  the  PKP  nodules 
are  less  than  those  of  the  other  nodule  types  because  their  potassium  oxide  content  of  0.86% 
stabilises  phlogopite  and  the  amphibole  k-richterite  (Table  Cl).  Both  minerals  contain  AI2O3, 
hence  reducing  the  proportion  of  garnet.  Furthermore  the  seismic  velocities  of  amphibole  and 
phlogopite  are  considerably  lower  than  are  those  of  garnet  (Table  Cl).  The  agreement  between 
the  PKP  parameters  and  those  of  the  best  fitting  seismic  model  therefore  suggests  that  mica, 
rather  than  garnet,  peridotite  is  the  rock  type  that  is  present  in  the  upper  mantle  over  much  of 
southern  Africa.  This  conclusion  is  consistent  with  the  petrological  observations.  Cox  (1983)  in 
particular  has  argued  that  the  potassium  concentration  in  the  Karoo  picrites  from  the  Lebombo 
monocline  requires  them  to  have  been  derived  from  two  end  members,  one  of  which  is  a  typical 
low  potassium  plume  basalt  and  the  other  a  K-rich  melt  with  4%  K2O  that  is  easily  produced 
by  melting  material  with  the  composition  of  PKP  nodules.  However,  as  Cox  (1983)  remarks,  the 


enormous  volume  of  the  picrites  require  a  large  source  of  the  K-rich  magma.  If  the  lid  consists  of 
PKP  material  and  is  60  km  thick,  1%  melting  generates  600  m  of  magma.  Cox’s  (1983)  model 
requires  the  2  km  thickness  of  the  picrites  to  contain  about  500  m  of  the  K-rich  melt.  There  is 
therefore  no  obvious  difficulty  in  accounting  for  the  geological  observations  in  this  way. 

The  difference  between  the  velocities  estimated  from  the  seismological  observations  below  a 
depth  of  122  km  and  those  from  the  deeper  nodules  is  considerably  larger  than  the  estimated 
uncertainties  (Fig.  10).  Though  presence  of  a  very  small  (  0.1%)  melt  fraction  with  a  dihedral 
angle  that  is  small  or  zero  seems  a  plausible  explanation  for  the  observations,  and  its  upward 
percolation  can  account  for  the  potassium  enrichment  seen  in  the  PKP  nodules,  it  is  not  likely 
that  the  whole  region  between  120  and  400  km  contains  melt.  Melting  is  generally  accompanied 
by  an  increase  in  volume,  and  hence  the  solidus  temperature  increases  with  pressure.  However, 
the  velocity  structure  of  the  LVZ  is  not  yet  well  constrained  by  the  seismic  observations  (see 
Fig.  10),  and  therefore  may  be  less  different  from  that  calculated  from  the  deeper  nodules  than 
our  model  suggests.  Furthermore,  all  the  nodules  consist  of  mantle  that  has  been  depleted  by 
considerable  amounts  of  basalt  extraction,  a  process  that  lowers  the  density  and  increases  the 
seismic  velocities.  The  upper  mantle  beneath  the  lithosphere  has  not  been  depleted  in  this  way, 
and  the  density  and  seismic  velocities  at  depths  greater  than  about  200  km  should  be  compared 
with  those  calculated  for  the  undepleted  material,  and  not  with  the  values  from  the  nodules. 

4.3  Comparison  with  Previous  Southern  African  Seismological  Models 

Figure  11  compares  our  earth  model  with  several  previously  published  models  for  southern 
Africa.  Bloch  et  al.  (1969)  used  fundamental  mode  Rayleigh  wave  phase  velocity  and  fundamental 
and  higher  mode  Love  and  Rayleigh  wave  group  velocities  to  determine  the  crust  and  upper  mantle 
shear  wave  structure  in  the  same  region  as  our  study,  and  we  used  their  mantle  model  for  our 
starting  model  in  the  waveform  inversion.  The  two  shear  wave  models  have  similar  upper  mantle 
lids  and  LVZ;  however,  the  average  mantle  lid  velocity  of  the  Bloch  et  al.  (1969)  model  is  higher 
than  in  our  model  and  the  average  velocity  in  the  LVZ  is  lower.  Our  lid  shear  wave  velocity  is 
controlled  by  the  Vs„  observations  (Fig.  5);  if  we  had  assumed  a  higher  Vs„,  this  would  result 
in  lower  velocities  in  the  LVZ.  Figure  12  compares  Bloch  et  al.  (1969)  dispersion  data  with  the 
dispersion  curves  computed  for  our  model.  The  observed  Rayleigh  wave  group  velocity  curves 
are  well  fit  by  the  group  velocity  curves  predicted  by  our  model  at  periods  greater  than  about  5 
seconds.  The  Love  wave  group  velocity  curves  are  well  fit  above  about  9  seconds  period  where 
the  fundamental  and  first  higher  mode  interfere.  The  predicted  Rayleigh  wave  phase  velocity 
curves  underestimate  the  observed  phase  velocities.  The  curvature  of  the  fundamental  mode 
Rayleigh  wave  phase  velocity  curves  near  40  s  period  (the  Airy  phase)  which  is  matched  in  Figure 
12,  provides  the  greatest  dispersion  constraint  on  the  lid  thickness  (Priestley  et  al,  1980).  The 
observed  phase  velocities  in  this  period  range  suggest  that  velocities  in  the  lower  part  of  the  lid 
may  be  slightly  higher  than  in  our  model;  however  Figure  7a  shows  that  synthetic  seismograms 
computed  for  the  Bloch  et  al.  (1969)  model  do  not  fit  the  observed  seismograms.  We  feel  that 
the  inversion  of  the  multi-mode  waveform  data  coupled  with  the  Vs„  data  from  Figure  5  provides 
stronger  constraints  on  the  upper  mantle  structure  than  is  provided  by  the  two-stage  process 
of  computing  the  fundamental  mode  dispersion  curve  and  inverting  that  for  the  upper  mantle 
structure. 

Cichowicz  <&£  Green  (1992)  use  the  partitioned  waveform  inversion  method  of  Nolet  (1990) 
to  invert  Rayleigh  wave  data  from  two  teleseismic  earthquakes  recorded  on  four  broadband  seis¬ 
mographs  deployed  along  a  profile  extending  southwest  from  the  central  Kaapvaal  Craton  into 
the  adjacent  Namaqua  mobile  belt.  Three  of  their  upper  mantle  shear  wave  models,  obtained 
by  including  3,  6,  and  11  eigenvectors  in  their  inversion,  are  shown  in  comparison  to  our  shear 
velocity  model  in  Figure  11.  The  model  with  3  eigenvectors  is  almost  horizontally  homogeneous 


whereas  the  other  models  show  an  upper  mantle  lid.  However  their  lid  is  significantly  thinner  and 
their  upper  mantle  shear  wave  LVZ  significantly  shallower  than  that  found  in  our  study.  Figure  13 
compares  reflectivity  synthetics  for  Event  2  computed  from  the  final  crustal  model  for  this  event, 
overlying  the  upper  mantle  models  of  Cichowicz  &  Green  (1992).  These  synthetic  waveforms 
show  little  resemblance  to  the  observed  waveforms.  We  used  the  Bloch  et  al.  (1969)  shear  wave 
structure  as  the  starting  mantle  model  in  our  inversion.  We  constructed  a  new  starting  model 
using  the  average  crustal  model  from  the  seismic  refraction  and  receiver  function  modeling,  and 
the  Cichowicz  Green  (1992)  n  =  6  mantle  model  and  reinverted  the  Event  2  seismograms  from 
SLR.  The  resulting  inversion  velocity  structure  which  employed  this  new  starting  model  converged 
towards  our  original  result  for  this  event,  giving  us  confidence  that  the  differences  between  our 
upper  mantle  model  and  those  of  Cichowicz  S/.  Green  (1992)  were  not  simply  the  result  of  our 
choice  of  a  starting  model. 

Clouser  &  Langston  (1990)  derived  upper  mantle  compressional  wave  models  for  both  cratonic 
and  mobile  belt  regions  of  southern  Africa  by  forward  modeling  waveforms  from  two  moderate 
earthquakes  in  Zimbabwe  recorded  on  long  period  seismographs  in  southern  Africa.  In  our  study 
of  the  Pal  waveforms  we  have  followed  the  same  approach  in  deriving  our  upper  mantle  P-wave 
velocity  model.  Clouser  &.  Langston  (1990)  conclude  that  a  model  with  an  upper  mantle  P-wave 
velocity  gradient  of  0.0033  s“^  and  a  P-wave  LVZ  no  shallower  than  170  km  best  fit  the  data 
for  the  craton.  This  velocity  model,  SACMOGc,  is  shown  in  Fig.  11  for  comparison  with  our 
model.  Since  SACMOGc  only  extends  to  275  km  depth,  we  have  extended  the  gradient  below  this 
depth.  The  seismograms  which  control  the  deeper  part  of  model  SACMOGc  were  recorded  at  an 
epicentral  distance  of  about  1100  km.  Nearer  events  could  be  fit  with  structures  containing  a 
P-wave  LVZ  as  shallow  as  130  km.  The  events  discussed  in  Appendix  B  have  epicentral  distances 
as  great  as  3000  km  and  hence  have  deeper  turning  depths.  Figure  14a  shows  that  reflectivity 
synthetics  computed  for  the  model  SACMOGc  closely  match  the  observed  Event  2  waveforms  at 
SLR  (A  =  1038  km).  However,  the  synthetics  P-wave  computed  for  SACMOGc  arrives  about 
3  s  earlier  than  the  observed  P-wave  arrival  time  of  the  March  9,  1989  earthquake  (Event  4, 
Table  1)  recorded  at  SLR  (A  =  1479  km)  (Fig.  14b).  This  early  arrival  time  arises  because  of 
the  steeper  upper  mantle  gradient  in  model  SACMOGc  compared  to  the  upper  mantle  gradient  in 
our  model.  However  the  narrow  P-wave  LVZ  in  model  SACMOGc  has  little  effect  on  the  deeper 
turning  P-wave  from  the  July  24,  1991  earthquake  (Event  10,  Table  1)  recorded  at  SUR  (A 
=  2106  km).  From  the  tests  we  have  made  in  comparing  synthetics  with  the  observations  (see 
Appendix  B),  we  conclude  that  if  a  P-wave  LVZ  exists  beneath  the  Archaean  core  it  is  deeper 
than  250  km  depth.  It  is  more  likely  that  there  is  no  P-wave  LVZ  beneath  the  Archaean  core. 

Figure  11  also  shows  the  P-wave  model  of  Green  (1978)  for  eastern  and  southern  Africa 
which  was  derived  from  travel  times  of  accurately  located  earthquakes  in  east  Africa  recorded  at 
seismographs  in  eastern  and  southern  Africa.  The  shallow  part  of  the  model  (<  250  km  depth)  is 
constrained  using  seismograms  of  east  African  earthquakes  recorded  at  east  African  stations  and 
therefore  characteristic  of  the  upper  mantle  beneath  the  East  African  Rift.  The  P-wave  velocity 
at  these  depths  is  low  compared  to  that  beneath  the  craton.  The  model  below  250  km  depth 
is  constrained  by  seismograms  recorded  on  the  shield  and  at  these  depths  Green’s  model  is  not 
significantly  different  from  the  model  we  have  determined  from  the  body  wave  modeling. 

4.4  Comparison  with  other  Earth  Models  for  Shield  Regions 

The  agreement  of  in  situ  velocities  and  density  estimates  at  upper  mantle  depths  beneath  the 
Kaapvaal  Craton  made  from  petrological  analysis  of  the  kimberlite  nodules  and  the  seismically 
determined  velocities  and  density  suggests  that  the  nodule  estimates  are  representative  of  the 
wider  region  of  the  Archaean  core.  It  is  interesting  to  compare  the  southern  African  results  with 
other  shield  regions  from  which  little  upper  mantle  nodule  data  is  available.  Figure  15  compares 


earth  models  for  the  Australian,  Baltic,  and  Canadian  shields  with  the  southern  African  model. 
The  shear  velocity  models  were  derived  from  either  forward  fitting  (Canadian  Shield  -  Brune  Si 
Dorman,  1963)  or  inversion  (Baltic  Shields  -  Calcagnile,  1991)  of  surface  wave  dispersion  data,  or 
by  forward  modeling  S-waves  recorded  in  the  regional  distance  range  (Kennett  et  al,  1994).  The 
compressional  velocity  models  were  determined  from  either  forward  modeling  long-period  (Baltic 
-  Given  Sl  Helmberger,  1980;  Canadian  Shields  -  LeFevre  Sl  Helmberger,  1989)  or  broadband 
(Australian  Shield  -  Bowman  Sl  Kennett,  1990)  P-waveforms  but  generally  at  larger  distance 
ranges  than  those  modeled  for  southern  Africa. 

The  upper  mantle  model  for  southern  Africa  has  features  similar  to  the  upper  mantle  beneath 
the  Canadian  and  Baltic  Shields.  The  shear  wave  structure  for  all  three  regions  shows  an  upper 
mantle  lid  with  shear  wave  velocities  varying  from  about  4.65  to  4.70  kms“T  Neither  the 
Canadian  nor  Baltic  shield  models  show  a  gradient  in  the  lid  but  this  is  probably  due  to  limitations 
in  the  method  of  analysis.  The  seismically  determined  velocities  in  the  lid  of  all  three  models 
are  in  reasonable  agreement  with  the  velocity  estimates  from  the  nodule  data  from  the  Kaapvaal 
Craton  kimberlites.  The  lid  overlies  a  shear  wave  LVZ  which  begins  at  about  110-120  km  depth. 
The  lowest  velocities  in  the  LVZ  are  at  depths  of  about  250  km  depth  in  each  of  the  three  regions. 
The  upper  part  of  the  Australian  shield  shear  wave  model  is  similar  to  the  shear  wave  structure 
for  the  other  three  shields.  However  the  S-wave  velocity  reversal  seen  below  the  other  shields  at 
120  to  150  km  depth  does  not  occur  beneath  the  Australian  shield  until  depths  of  200  to  225 
km. 

There  is  more  variation  in  the  P-wave  models  for  the  four  shield  regions.  Below  about  220  km 
depth  the  Canadian  and  Baltic  Shield  models  are  similar  to  the  southern  African  model;  however, 
above  this  the  models  are  significantly  different.  Both  the  Canadian  and  Baltic  Shield  models 
have  a  high  velocity  P-wave  lid  above  a  P-wave  LVZ  but  the  nature  of  the  lid  and  the  LVZ  is 
different  in  the  two  models.  The  Pn  velocity  in  the  Baltic  Shield  model  is  about  8.2  kms~^  and 
the  lid  extends  to  a  depth  of  about  140  km;  the  Canadian  Shield  model  has  a  Pn  velocity  of 
about  8.4  kms“^  and  the  lid  extends  to  a  depth  of  about  170  km.  In  contrast,  the  southern 
African  Pn  velocity  is  8.1  kms”^  and  there  is  no  P-wave  LVZ.  The  P-wave  velocity  gradients 
in  the  lid  of  the  Canadian  and  Baltic  Shield  models  are  similar  to  that  in  the  southern  African 
model.  The  P-wave  velocities  of  the  Canadian  and  Baltic  Shield  models  do  not  agree  with  the 
estimates  from  the  kimberlite  nodules  from  the  Kaapvaal  craton.  The  Australian  Shield  P-wave 
model  does  not  contain  a  LVZ  but  has  a  significantly  higher  velocity  than  the  southern  African 
model  and  does  not  correspond  to  the  same  structure  as  the  other  models  below  220  km  depth. 

5.  CONCLUSIONS 

We  present  a  new  1— D  velocity  and  density  model  for  the  crust  and  upper  mantle  beneath 
southern  Africa.  The  crustal  velocity  structure  is  derived  from  modeling  teleseismic  receiver 
function  data  and  from  published  seismic  refraction  crustal  models.  The  mantle  model  is  derived 
from  inversion  of  multi-mode  regional  surface  waveforms  and  forward  modeling  of  Pn/  waveforms. 
The  average  crust  and  upper  mantle  model  derived  from  analysis  of  data  from  twelve  paths  largely 
confined  to  the  Archaean  core  of  southern  Africa  has  the  following  major  features.  The  42  km 
thick  crust  consists  of  four  layers:  a  2  km  thick  5.05  kms~^  surface  layer,  a  5  km  thick  6.08 
kms"^  upper  crust,  a  20  km  thick  6.30  kms-^  mid-crust,  and  a  15  km  thick  6.73  kms"^  lower 
crust.  Below  the  crust  there  is  an  80  km  thick  upper  mantle  lid.  Vp^  and  Vs„  are  8.09  and  4.62 
kms“^,  respectively,  and  the  compressional  and  shear  velocity  gradients  in  the  lid  are  0.0008  s“^ 
and  0.0013  s■^  respectively.  Below  the  lid  there  is  a  substantial  S-wave  LVZ  with  the  shear 
wave  velocity  dropping  to  4.32  kms”^  at  250  km  depth.  The  S-wave  LVZ  is  required  to  both 
explain  the  Vs^  measurements  and  fit  the  regional  surface  waveforms.  Forward  modeling  of  Pn/ 
waveforms  does  not  require  an  upper  mantle  P-wave  LVZ.  Below  125  km  depth  the  P-wave 
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gradient  increases  to  0.0015  s  ^  and  increases  again  to  0.0035  between  250  km  depth  and 
the  410  discontinuity.  The  distribution  of  earthquakes  and  seismographs  used  in  this  study  has 
not  provided  adequate  data  to  resolve  the  question  of  deep  lithospheric  differences  in  the  structure 
beneath  the  Archaean  cratons  and  the  Proterozoic  mobile  belts. 

We  have  estimated  the  upper  mantle  lid  velocities  and  density  from  the  composition  of  the 
nodules  brought  up  by  kimberlites  through  the  Kaapvaal  Craton.  These  values  are  in  good 
agreement  with  the  seismically  determined  upper  mantle  lid  velocities  and  density  for  the  Archaean 
core  of  southern  Africa,  suggesting  that  the  nodule  data  is  representative  of  the  wider  region  and 
not  just  the  region  beneath  where  they  are  found. 
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Table  1.  Regional  Events  in  Southern  Africa 


No 

DATE 

TIME 

LAT. 

LONG. 

DEPTH 

Mb 

Ms 

H 

1983 

7 

7 

20 

35 

37.4 

-7.41 

27.88 

10 

5.8 

5.5 

El 

1986 

7 

18 

15 

7 

52 

-16.36 

28.48 

17 

5.4 

4.8 

1988 

4 

16 

21 

17 

9.9 

-10.26 

27.71 

10 

5.4 

4.7 

4+ 

1989 

3 

9 

2 

37 

0.2 

-13.74 

34.41 

30 

5.6 

5.4 

5+ 

1989 

3 

10 

21 

49 

45.8 

-13.70 

34.45  . 

30 

6.0 

6.1 

6+ 

1989 

9 

29 

6 

16 

31 

-30.60 

28.70 

7 

4.8 

4.1 

V 

1990 

9 

26 

23 

8 

23.1 

-28.10 

26.77 

2 

5.3 

4.3 

8+ 

1990 

11 

3 

0 

31 

28.9 

-21.35 

33.28 

10 

4.8 

5.1 

9t 

1991 

5 

10 

1 

12 

38.5 

-17.35 

24.98 

10 

4.8 

3.8 

10+ 

1991 

7 

24 

13 

54 

51.9 

-18.27 

34.86 

32 

5.1 

4.9 

11+ 

1991 

11 

3 

18 

28 

55 

-26.89 

26.64 

1 

5.0 

12+ 

1994 

8 

18 

0 

45 

47.2 

-7.43 

31.75 

25 

6.0 

5.7 

13+ 

1994 

10 

30 

6 

6 

27.4 

-28.03 

26.74 

5 

5.6 

4.7 

14+ 

1995 

1 

2 

15 

16 

2.2 

-26.36 

27.22 

5 

4.3 

t|SC 

*PDE 


Table  2.  Regional  Earthquake  Source  Mechanisms 


CMT  solutions 

k 

source  mechanisms  used 

Event 

depth 

Mo -1024 

NPl 

(km) 

(dyn  cm) 

6 

A 

(t> 

8 

A 

830707 

24.1 

5.700 

200 

41 

-111 

6.0 

5.700 

200 

41 

-111 

860718 

15.0 

0.783 

238 

43 

-84 

12.0 

0.783 

238 

43 

-84 

880416 

15.0 

0.535 

10 

26 

-120 

10.0 

0.428 

353 

18 

-139 

890309 

15.0 

3.690 

128 

23 

-118 

30.0 

1.845 

128 

23 

-118 

890310 

15.0 

31.100 

142 

17 

-102 

30.0 

15.552 

154 

34 

-92 

901103 

15.0 

0.592 

345 

35 

-111 

10.0 

0.700 

310 

20 

-132 

910724 

24.7 

0.518 

180 

45 

-90 

32.0 

0.518 

180 

45 

-90 

940818 

30.9 

9.300 

306 

28 

-110 

30.9 

9.300 

306 

28 

-no 

*Harvard  centroid  moment  tensor  catalogue 
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Table  3.  Southern  African  Digital  Seismograph  Stations. 


Code 

Lat. 

Long. 

Elev. 

Site- Name 

SLR 

-25.7349 

28.2816 

1348.0 

Silverton,  South  Africa 

SUR 

-32.3797 

20.8117 

1770.0 

Sutherland,  South  Africa 

BOSA 

-28.6131 

25.4156 

1280.0 

Boshof,  South  Africa 

LBTB 

-25.0145 

25.5970 

1128.0 

Lobatse,  Botswana,  Africa 

LSZ 

-15.2766 

28.1882 

1184.8 

Lusaka,  Zambia,  Africa 

TSUM 

-19.2022 

17.5838 

1240.0 

Tsumeb,  Namibia,  Africa 

Table  4.  Earth  Model  for  July  18,  1986  -  SLR  Path 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

kms“^ 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

kms~^ 

Vs 

km 

P 

gcm“^ 

0 

1 

4.52 

2.44 

2.09 

199 

21 

8.34 

4.26 

3.41 

1 

5 

5.80 

3.67 

2.31 

219 

21 

8.38 

4.14 

3.44 

6 

20 

6.45 

3.71 

2.80 

240 

23 

8.44 

4.10 

3.45 

26 

14 

6.60 

3.92 

2.88 

263 

23 

8.53 

4.10 

3.46 

40 

13 

8.10 

4.64 

3.35 

287 

23 

8.62 

4.10 

3.48 

53 

13 

8.12 

4.65 

3.35 

310 

17 

8.69 

4.17 

3.49 

66 

13 

8.13 

4.67 

3.35 

327 

17 

8.75 

4.29 

3.51 

78 

13 

8.15 

4.68 

3.35 

343 

17 

8.81 

4.42 

3.52 

91 

13 

8.17 

4.69 

3.35 

360 

17 

8.88 

4.55 

3.54 

104 

13 

8.20 

4.72 

3.35 

377 

17 

8.94 

4.68 

3.55 

116 

21 

8.21 

4.56 

3.35 

393 

17 

9.03 

4.87 

3.57 

137 

21 

8.24 

4.48 

3.37 

410 

42 

9.36 

5.07 

3.76 

158 

21 

8.27 

4.41 

3.38 

452 

8 

9.50 

5.16 

3.81 

178 

21 

8.31 

4.33 

3.40 

50 


Table  5.  Average  Velocity  Model  for  Southern  Africa 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

kms~^ 

P 

g  cm~^ 

Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

kms“^ 

P 

g  cm“^ 

0 

2 

5.05 

2.88 

2.36 

232 

10 

8.33 

4.32 

3.43 

2 

5 

6.08 

3.59 

2.74 

242 

10 

8.35 

4.32 

3.43 

7 

20 

6.30 

3.68 

2.88 

252 

10 

8.38 

4.32 

3.44 

27 

15 

6.73 

4.04 

3.00 

262 

10 

8.41 

4.34 

3.45 

42 

10 

8.09 

4.62 

3.30 

272 

10 

8.44 

4.35 

3.45 

52 

10 

8.09 

4.63 

3.31 

282 

10 

8.48 

4.38 

3.46 

62 

10 

8.10 

4.65 

3.32 

292 

10 

8.51 

4.41 

3A7 

72 

10 

8.11 

4.66 

3.32 

302 

10 

8.55 

4.43 

3.48 

82 

10 

8.12 

4.68 

3.33 

312 

10 

8.58 

4.47 

3.48 

92 

10 

8.13 

4.70 

3.34 

322 

10 

8.62 

4.50 

3.49 

102 

10 

8.14 

4.71 

3.35 

332 

10 

8.65 

4.54 

3.50 

112 

10 

8.15 

4.72 

3.35 

342 

10 

8.69 

4.58 

3.51 

122 

10 

8.17 

4.50 

3.37 

352 

10 

8.73 

4.61 

3.51 

132 

10 

8.18 

4.46 

3.37 

362 

10 

8.76 

4.65 

3.52 

142 

10 

8.20 

4.44 

3.38 

372 

10 

8.80 

4.68 

3.53 

152 

10 

8.21 

4.43 

3.38 

382 

10 

8.83 

4.72 

3.54 

162 

10 

8.22 

4.41 

3.39 

392 

10 

8.87 

4.76 

3.54 

172 

10 

8.24 

4.39 

3.39 

402 

9 

8.90 

4.79 

3.55 

182 

10 

8.25 

4.38 

3.40 

411 

10 

9.31 

5.01 

3.76 

192 

10 

8.26 

4.36 

3.40 

421 

10 

9.34 

5.04 

3.78 

202 

10 

8.28 

4.35 

3.41 

431 

10 

9.38 

5.06 

3.79 

212 

10 

8.30 

4.34 

3.42 

441 

10 

9.42 

5.09 

222 

10 

8.31 

4.33 

3.42 

451 

0 

9.45 

5.11 

3.81 

Table  6.  Upper  Mantle  Parameters  for  Nodule  Analysis. 


Nodule 

Comp. 

Temp. 

Press. 

Depth 

Vp 

Vs 

P 

PKP 

Sp.  P. 

371 

1.61865 

50 

8.078 

4.627 

3.333 

PKP 

Sp.  P. 

557 

2.42798 

75 

8.086 

4.632 

3.338 

PKP 

Gar.  P 

743 

3.23730 

100 

8.094 

4.637 

3.343 

PKP 

Gar.  P 

854 

3.72290 

115 

8.099 

4.641 

3.346 

Depleted 

Gar.  P 

1058 

4.04663 

125 

8.327 

4.830 

3.298 

Fertile 

Gar.  P 

1380 

5.82714 

180 

8.394 

4.864 

3.311 

51 


FIGURE  CAPTIONS 


Figure  1.  Tectonic  map  of  southern  Africa  showing  the  major  geologic  subdivisions  (after  Cahen 
&  Snelling,  1966).  The  locations  of  the  kimberlite  pipes  are  shown  by  the  solid  diamonds  and 
are  those  from  which  the  samples  used  by  Tainton  ii  McKenzie  (1994)  were  obtained. 

Figure  2.  Source  -  receiver  paths  for  regional  earthquake  seismograms  used  in  this  study  superim¬ 
posed  on  the  major  crustal  subdivisions  of  southern  Africa.  Solid  triangles  denote  seismograph 
stations;  stars  denote  earthquake  locations;  solid  lines  denote  paths  for  the  events  used  in  the 
regional  waveform  modeling;  and  dotted  lines  denote  paths  for  events  used  in  constructing 
the  Pn,  Sn,  Lg  composite  record-section.  Events  2-SLR,  5-SLR,  8-SLR,  and  12-LSZ  were 
used  for  waveform  modeling  and  the  record  section.  The  numbers  for  the  earthquakes  refer 
to  Table  1. 

Figure  3.  Results  of  the  teleseismic  receiver  function  analysis  for  data  from  the  station  BOSA. 
The  solid  line  is  the  observed  receiver  function  and  the  dashed  line  denotes  ±  one  standard 
deviation.  The  dotted  line  corresponds  to  the  synthetic  receiver  function  computed  from  the 
crustal  velocity  structure  shown  in  the  inset.  We  have  annotated  the  receiver  function  with 
the  phases  corresponding  to  the  relevant  features  of  the  model  shown  in  the  inset.  In  labeling 
these  phases  we  have  used  the  convention  of  Bath  &.  Stefansson  (1966).  The  I  and  m  in  the 
inset  velocity  model  denote  the  boundary  where  the  corresponding  P-to-S  conversions  in  the 
synthetic  receiver  function  occur. 

Figure  4.  Measured  values  of  crustal  thickness  in  southern  Africa.  Large  circles  denote  location  of 
teleseismic  receiver  function  crustal  thickness  estimates  from  this  study;  small  circles  denote 
location  of  seismic  refraction  crustal  thickness  estimates.  The  letter  in  parenthesis  following 
the  refraction  estimates  corresponds  to  the  following  references:  (i)  Willmore  et  al,  1952,  (ii) 
Gane  et  al,  1956,  (iii)  Hales  &  Sacks,  1959,  (iv)  Baier  et  al,  1983,  (v)  Stuart  ii  Zengeni, 
1987,  (vi)  Green  ii  Durrheim,  1990,  (vii)  Durrheim  ii  Green,  1992,  (viii)  Durrheim  et  al, 
1992,  (ix)  Wright  ii  Hall,  1990. 

Figure  5.  High-frequency  body  wave  composite,  reduced  record-section  from  seismograms  of  the 
earthquakes  and  seismographs  shown  in  Figure  2.  The  travel  time  relationships  for  the  P^, 
Sn,  and  Lg  phases  are  shown  on  the  right  hand  edge  of  the  plot. 

Figure  6.  Vertical  component  short-period  and  three  component  long-period  seismograms  for 
the  July  18,  1986  (Event  2,  Table  1)  earthquake  recorded  at  SLR.  The  seismograms  are 
annotated  with  the  predicted  lASPEl  P  and  S  wave  arrival  times  and  with  our  interpretation 
of  the  various  arrivals. 

Figure  7.  (a)  Comparison  of  the  July  18,  1986  earthquake  (Event  2,  Table  1)  long-period 
seismograms  recorded  at  SLR  (solid  line)  with  normal  mode  synthetics  seismograms  computed 
for  the  starting  model,  the  CMT  source  mechanism,  and  the  ISC  depth  (dotted  line),  (b) 
Comparison  of  the  same  seismograms  and  locked-mode  synthetic  seismograms  (dotted  lines) 
computed  using  the  velocity  model  obtained  by  inversion.  For  this  comparison  the  CMT 
source  mechanism  was  also  used  but  the  focal  depth  was  reduced  to  12  km.  (c)  Comparison 
of  the  starting  S-wave  velocity  model  (dashed  line)  and  the  S-wave  inversion  model  (solid 
line).  Qualitative  bounds  on  our  inversion  model  derived  from  the  test  discussed  in  Appendix 
B  are  shown  as  dotted  lines. 

Figure  8.  Vertical  component  short-period  and  three  component  long-period  P-wave  seismo¬ 
grams  for  Event  2.  52 


Figure  9.  (a)  Comparison  of  long-period  P-waveform  of  Event  2  (solid  line)  with  reflectivity 
synthetic  seismograms  computed  for  the  inversion  (dashed  line),  and  for  the  final  (dotted 
line)  velocity  models.  The  Pi  waveform  which  is  primarily  controlled  by  the  crustal  structure 
is  equally  well  fit  by  synthetics  waveforms  for  both  models.  However,  the  amplitude  of  the 
first  arrival  is  better  fit  by  the  synthetic  computed  for  the  final  model,  (b)  Comparison  of 
long  period  seismograms  for  Event  2  with  whole  wavefield  reflectivity  synthetic  seismograms 
for  the  final  earth  model,  (c)  Earth  model  from  surface  wave  inversion  (dotted  lines)  and  final 
earth  model  (solid  line)  from  forward  modeling  waveforms  for  this  Event  2.  The  shear 
wave  velocity  model  remained  fixed  during  this  stage  of  the  analysis.  The  P-  and  S-wave 
velocity  is  in  kms”^  and  the  density  is  in  gmc"^. 

Figure  10.  Average  earth  model  (solid  lines)  for  southern  Africa  determined  by  computing  the 
mean  and  standard  deviation  of  the  velocities  and  density  at  each  depth  for  the  twelve  path- 
average  earth  models.  The  dotted  lines  denote  the  range  of  velocity  models  from  analysis  of 
the  individual  source— receiver  paths.  The  solid  circles  denote  velocity  and  density  estimates 
from  geochemical  analysis  of  the  garnet  peridotite  nodules  from  kimberlites  on  the  Kaapvaal 
Craton.  The  x-axis  scale  is  in  kms“^  for  the  P-  and  S-wave  velocities  and  gmcm“^  for  the 
density. 

Figure  11.  Southern  African  velocity  model  from  this  study  (solid  line)  compared  with  previous 
models  for  southern  Africa:  (a)  shear  wave  model  of  Bloch  et  al.  (1969)  (dotted  line), 
compressional  wave  model  of  Clouser  Langston  (1990)  (dotted  line),  compressional  wave 
model  of  Green  (1978)  (short  dashed  line);  (b)  shear  wave  models  of  Cichowicz  8i  Green 
(1992)  (short  dashed  line,  n=3;  long  dashed  lines,  n=6;  dash-dot  line,  n=ll). 

Figure  12.  Comparison  of  the  fundamental  and  higher  mode  Love  (a)  and  Rayleigh  (b)  wave 
group  velocity  data  and  fundamental  mode  Rayleigh  (b)  wave  phase  velocity  data  from  Bloch 
et  al.  (1969)  (dotted  line),  and  dispersion  curves  (solid  lines)  computed  for  the  velocity  and 
density  model  of  Figure  10  and  Table  5. 

Figure  13.  Comparison  of  Event  2  long-period  seismograms  (solid  line)  with  reflectivity  synthetics 
computed  for  the  upper  mantle  shear  wave  models  for  the  Kaapvaal  Craton  from  Cichowicz 
&  Green  (1992)  n=6  (dotted).  This  fit  should  be  compared  with  the  fit  of  the  reflectivity 
synthetic  seismograms  computed  for  the  final  earth  model  (Table  4)  for  this  event  shown  in 
Figure  9. 

Figure  14.  (a)  Comparison  of  observed  vertical  component  waveforms  (solid  lines)  for  the  July 
18,  1986  earthquake  (Event  2,  Table  1)  recorded  at  SLR  with  synthetic  waveforms  computed 
for  the  final  earth  model  for  this  event  (dashed  line)  and  for  the  SACM06c  P-wave  model 
of  Clouser  8i  Langston  (1990)  (dotted  line),  (b)  Same  comparison  as  in  (a)  but  for  the 
seismogram  for  the  March  9,  1989  earthquake  (Event  4,  Table  1)  recorded  at  SLR.  (c)  Same 
comparison  as  in  (a)  but  for  the  seismograms  of  the  July  24,  1991  earthquake  (Event  10, 
Table  1)  recorded  at  SDR. 

Figure  15.  Southern  African  velocity  model  from  this  study  (solid  line)  compared  with  veloc¬ 
ity  models  for  other  cratonic  regions:  (i)  Canadian  shield  (long  dashed  line)  (LeFevre  & 
Helmberger,  1989;  Brune  Dorman,  1963),  (ii)  Australian  shield  (dotted  line)  (Bowman  & 
Kennett,  1990;  Kennett  et  al.,  1994),  and  (iii)  the  Baltic  shield  (short  dashed  line)  (Given 
Helmberger,  1980;  Calcagnile,  1991). 
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FIGURE  7b 
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APPENDIX  A  Teleseismic  Receiver  Function  Analysis 

Seismograms  of  teleseismic  P-waveforms  contain  the  combined  effects  of  the  earthquake 
source,  propagation  through  the  Earth,  and  the  instrument  response.  To  use  these  seismic  phases 
to  determine  earth  structure  beneath  a  broadband  seismograph,  the  effects  of  the  crust  and  upper 
mantle,  and  the  site  response  or  receiver  function,  must  be  isolated  from  these  other  factors.  We 
have  computed  radial  and  tangential  receiver  functions  by  computing  the  spectral  ratio  of  the 
radial  or  tangential  component  spectrum  and  the  vertical  component  spectrum  (Langston,  1979; 
Ammon,  1991).  The  deconvolution  is  stabilized  by  specifying  a  minimum  spectral  trough  fill  c 
for  the  vertical  spectrum  to  avoid  noise  generated  by  dividing  by  small  values.  The  resulting 
receiver  function  is  smoothed  by  a  Gaussian,  a,  which  limits  the  high  frequency  content.  All 
deconvolutions  in  this  study  were  computed  for  c  =  0.001  and  a  =  1.0  which  were  found  to 
produce  the  most  stable  set  of  deconvolutions  for  the  whole  P-wave  data  set.  This  Gaussian 
limits  the  frequency  content  of  the  receiver  functions  to  frequencies  below  about  0.4  Hz,  resulting 
in  reduced  resolution  of  crustal  velocity  features.  However,  these  low  frequency  receiver  functions 
are  able  to  resolve  the  main  crustal  features  for  comparison  with  the  seismic  refraction  results 
where  they  exist,  and  provide  a  gross  crustal  structure  estimate  where  seismic  refraction  results 
do  not  exist.  The  spectral  ratio  is  then  transformed  back  into  the  time  domain  giving  a  time 
display  of  the  P-to-S  conversions  caused  by  the  velocity  impedance  contrast  beneath  the  site. 

Figure  A1  shows  single  event  and  stacked  receiver  functions  for  BOSA.  The  individual  radial 
receiver  functions  are  dominated  by  a  large  amplitude  arrival  4  s  after  direct  P.  There  are  two 
smaller  amplitude  arrivals  in  each  of  the  individual  radial  receiver  functions;  a  positive  arrival  ~ 
14  seconds  after  direct  P  and  a  negative  arrival  ~  18  seconds  after  direct  P.  We  have  stacked 
individual  receiver  functions  for  events  occurring  over  a  small  range  of  ray  parameter  and  azimuth 
and  computed  the  one-standard  deviation  bounds.  The  4,  14,  and  18  second  arrivals  in  the 
stacked  radial  receiver  function  are  large  compared  with  the  amplitude  of  the  stacked  tangential 
receiver  function  at  the  same  times.  However  there  are  coherent  tangential  arrivals  at  6  and  19 
seconds  which  may  possibly  interfere  with  these  arrivals.  If  we  interpret  the  4  second  phase  as 
a  P-to-S  conversion  at  the  Moho  (Ps)  and  assume  an  average  crustal  P-wave  velocity  of  6.4 
km  and  a  Poisson’s  ratio  of  0.25,  the  time  of  this  arrival  suggests  a  35  km  thick  crust  beneath 
BOSA. 

The  stacked  radial  receiver  function  was  first  forward  modeled  using  a  simple  single  layer 
crust  over  a  half-space  model  based  on  the  seismic  refraction  derived  P-wave  crustal  model  of 
Durrheim  &.  Green  (1992)  with  shear  wave  velocities  assuming  a  Poisson’s  ratio  of  0.25,  and 
densities  assuming  the  relationship  p  =  0.77  +  0.32Vp.  Once  the  major  features  of  the 
receiver  function  were  matched  by  this  simple  modeling,  the  stacked  radial  receiver  function  was 
inverted  using  a  time  domain,  linearized  inversion  procedure  (Ammon  et  al,  1990)  using  a  more 
complex  starting  earth  model  based  the  Durrheim  &.  Green  (1992)  model.  The  starting  model 
was  parameterized  as  a  stack  of  1  to  3  km  thick  layers  to  a  depth  of  70  km.  This  model  was 
randomly  perturbed  into  24  different  starting  models  and  the  radial  receiver  function  inverted  for 
a  1-D  velocity  model  of  the  receiver  structure  by  minimizing  the  difference  between  the  observed 
radial  receiver  function  and  a  synthetic  predicted  for  the  structure.  This  pseudo-Monte  Carlo 
approach  enabled  us  to  reduce  the  dependence  of  the  convergence  of  the  inversion  model  on 
the  form  of  the  starting  model.  Once  the  depths  of  the  significant  seismic  impedance  contrast 
were  identified,  the  initial  model  was  reparameterized  with  thicker  layers  and  the  radial  receiver 
function  re-inverted.  Using  this  approach  we  have  determined  a  crustal  model  beneath  each 
station  with  the  minimum  number  of  layers  which  give  a  satisfactory  fit  to  the  data.  Since  the 
receiver  function  consists  of  few  events,  we  concentrated  on  fitting  the  early  part  of  the  data  and 
did  not  emphasize  fitting  later  arrivals.  The  results  for  BOSA  are  shown  in  Figure  3  and  discussed 
in  the  main  text.  Discussion  of  the  receiver  function  analysis  for  the  other  sites  is  found  in  the 
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captions  to  Figures  A2-5.  The  format  for  these  figures  is  the  same  as  that  of  Figure  3. 

We  would  like  to  re-emphasize  that  the  purpose  of  analyzing  the  receiver  function  data  was 
not  to  produce  definitive  models  of  the  crustal  velocity  structure  beneath  the  four  stations  but 
to  provide  some  control  on  the  gross  crustal  structure  (average  velocity  and  Moho  depth)  to  use 
in  the  modeling  of  the  regional  waveforms. 

APPENDIX  A  -  FIGURE  CAPTIONS 

Figure  Al.  Teleseismic  receiver  functions  for  BOSA  from  four  deep  earthquakes  in  South  America 
(A  =  79°  ±  3°,  backazimuth  =  244°  ±  5°).  The  four  upper  panels  are  individual  receiver 
functions  and  the  two  bottom  panels  are  the  stacked  radial  and  transverse  receiver  functions. 
For  the  individual  receiver  functions  the  solid  line  denotes  the  radial  and  the  dotted  line  denotes 
the  transverse;  for  the  stack  receiver  function  the  solid  line  denotes  the  receiver  function  and 
the  dotted  line  denotes  the  ±  1  cr  bounds. 

Figure'A2.  This  figure  and  Figures  A3-5  are  all  in  the  same  format  as  Figure  3.  The  LBTB 
receiver  function  consists  of  a  stack  of  two  events  from  the  Hindu  Kush  region  (A  =  75°±1°, 
backazimuth  =  37°  ±  1°).  The  Psm  arrival  at  LBTB  follows  the  direct  arrival  by  almost  5 
seconds,  about  a  second  longer  than  the  Ps  -  P  time  at  BOSA.  Assuming  the  same  crustal 
velocities  in  this  part  of  the  craton,  the  LBTB  Ps  -  P  time  difference  implies  a  crust  r\j  8-10 
km  thicker  beneath  LBTB  than  beneath  BOSA.  However  there  is  a  significant  arrival  at  6 
seconds  on  the  tangential,  indicating  an  off-azimuth  arrival.  This  may  cause  distortion  in  the 
Ps  arrival.  There  is  a  second  negative  arrival  at  about  22  seconds  which  could  be  interpreted 
as  the  PpSms  +  PsPms  phase.  The  inversion  result  shows  a  relatively  high  shear  velocity  in 
the  surface  layer  and  a  nearly  uniform  3.77  kms“^  crust  between  4  and  34  km  depth.  The 
velocity  increases  to  3.95  kms"^  at  34  km  depth  and  there  is  a  sharp  crust  -  mantle  transition 
at  44  km  depth  where  the  velocity  increases  from  4.13  to  4.53  kms“^  The  synthetic  receiver 
function  for  this  velocity  model  fits  the  timing  of  the  main  features  of  the  observed  receiver 
function  but  the  amplitudes  are  not  well  fit.  The  coherency  is  poor  for  PpPfe  and  PpPms 
which  may  partially  result  from  interference  from  off-azimuth  energy  arriving  at  the  same 
time  as  indicated  by  the  large  coherent  arrival  on  the  tangential  component  between  12  and 
18  seconds.  The  amplitude  of  this  tangential  arrival  is  much  larger  than  the  synthetic  radial 
arrival  at  about  the  same  time.  The  timing  of  the  PpSms  +  PsPms  synthetic  and  observed 
arrival  agrees  but  there  is  a  large  discrepancy  in  amplitude.  While  the  correlation  between 
the  observed  and  synthetic  receiver  function  is  not  as  good  as  we  would  like,  the  Ps-P  time 
difference  does  suggest  that  the  crust  beneath  this  part  of  the  craton  is  thicker  than  to  the 
south  beneath  BOSA. 

Figure  A3.  The  SUR  receiver  function  stack  consists  of  three  events  from  South  America  (A  = 
75°  ±  6°,  backazimuth  =  253°  ±  6°).  The  Psm  arrival  follows  direct  P  by  nearly  6  seconds 
indicating  either  a  thicker  crust  or  lower  velocity  crust  compared  to  the  region  to  the  NE. 
There  are  no  other  clear  phases  in  the  radial  receiver  function.  The  velocity  model  deduced 
from  this  data  consists  of  a  nearly  uniform  crust  below  ~  5  km  depth,  with  a  gradual  velocity 
increase  starting  at  about  35  km  depth  to  mantle  velocities  below  51  km  depth.  The  synthetic 
receiver  function  for  this  model  provides  a  reasonably  good  fit  to  the  main  features  of  the 
observed  data. 

Figure  A4.  This  SUR  receiver  function  consists  of  a  stack  of  two  events  from  the  South  Sandwich 
Islands  (A  =  40°  ±  1°,  backazimuth  =  219°  ±  1°).  The  radial  receiver  function  for  this 
azimuth  shows  an  oscillatory  behavior  which  persists  for  10  seconds  following  Ps.  This  is  a 
common  feature  of  both  of  the  single  event  receiver  functions  comprising  the  stack  but  it  is 
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not  clear  how  these  oscillations  arise.  Except  for  the  coherent  arrival  at  8  seconds,  there  is 
little  indication  of  ofF-azimuth  arrivals  which  may  contaminate  the  radial  receiver  function. 
The  inversion  results  for  the  receiver  function  from  this  azimuth  are  similar  to  those  from  the 
South  American  azimuth  at  SUR  except  that  the  Moho  in  this  model  is  at  53  km  depth. 

Figure  A5.  The  receiver  function  for  LSZ  consists  of  a  single  event  from  South  America(A  =  85°, 
backazimuth  =  244°).  The  Psm  arrival  follows  the  direct  arrival  by  4  seconds.  Later  large 
arrivals  on  the  radial  component  are  questionable  because  of  the  large  amplitude  arrivals  seen 
on  the  tangential  component  starting  at  about  7  seconds.  The  inversion  result  shows  a  crust 
with  a  strong  gradient  in  the  shallow  crust,  a  uniform  3.63  kms”^  upper  crust  extending 
from  4  to  25  km  depth,  and  a  gradational  lower  crust  with  velocities  increasing  from  3.87 
kms-^  at  25  km  depth  to  4.2  kms"^  at  40  km  depth.  The  transition  from  lower  crust  to 
mantle  velocities  occurs  over  a  2  km  transition.  The  synthetic  Ps  phase  underestimates  the 
amplitude  of  the  observed  Ps  phase. 
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FIGURE  A4 
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APPENDIX  B  Regional  Waveform  Analysis 

This  appendix  discusses  in  more  detail  the  determination  of  path-average  velocity  and  density 
structure  from  the  regional  seismic  waveforms  of  the  eight  earthquakes  listed  in  Table  1.  We  give 
additional  details  of  the  analysis  procedure  we  have  followed  using  the  July  18,  1986  event  as  an 
example,  and  then  summarize  the  results  for  the  other  events. 

The  seismograms  of  the  July  18, 1986  earthquake  recorded  at  SLR  have  been  described  above. 
The  earthquake  occurred  within  the  Irumide  mobile  belt  north  of  the  Zimbabwe  Craton  and  the 
propagation  path  was  largely  confined  to  the  Archaean  core.  The  epicenter  was  almost  due  north 
of  SLR  so  the  recorded  seismograms  are  close  to  normally  rotated  (E-W  tangential,  N-S  radial). 
We  first  compare  the  observed  waveform  with  synthetic  waveforms  (Fig  7a)  computed  using 
the  CMT  mechanism  and  an  earth  model  consisting  of  the  average  crustal  model  determined 
above  (section  3.2)  overlying  the  mantle  shear  wave  velocity  structure  determined  by  Bloch  et 
al  (1969)  from  surface  wave  dispersion  analysis.  This  earth  model  is  parameterized  in  terms  of 
flat-lyfng  layers  of  constant  Vs;  Vp  is  related  to  Vs  by  assuming  a  Poisson’s  ratio  a  of  0.25  and 
the  density  is  related  to  Vp  by  />  =  0.77  +  0.32  Vp.  The  compressional  and  shear  wave  Q's  were 
assumed  to  be  1000  and  500,  respectively,  for  the  crust  (Clouser  &  Langston,  1990)  and  700  and 
300,  respectively,  for  the  mantle.  We  then  altered  the  source  mechanism  by  ±10°  in  strike,  dip, 
and  slip,  and  the  hypocenter  by  ±10  km  in  depth  and  distance  to  verify  that  the  observed  misfit 
was  not  the  result  of  errors  in  the  hypocenter  or  source  mechanism.  We  could  not  improve  the 
synthetic  fit  to  the  observed  waveform  in  this  way. 

To  invert  the  waveform  data  for  earth  structure  we  first  parameterized  the  starting  model 
in  terms  of  a  large  number  of  thin,  homogeneous  layers.  We  constrain  the  inversion  using 
bounds  on  the  crustal  structure  from  the  seismic  refraction  and  teleseismic  receiver  function 
studies  (Fig.  4)  and  bounds  on  the  velocity  of  the  upper  portions  of  the  mantle  lid  from  the 
Pn  and  Sn  traveltime  data  (Fig.  5).  We  then  inverted  the  three-components  seismograms  using 
the  differential  seismogram  inversion  technique  of  Gomberg  &  Masters  (1988).  Once  the  main 
character  of  the  model  boundaries  became  apparent  from  the  initial  inversions,  the  number  of 
layers  was  reduced  by  grouping  layers  with  similar  parameters  thus  fixing  the  main  boundaries; 
the  final  velocities  and  density  of  these  thicker  layer  units  were  determined  by  reinverting  the 
waveforms.  We  have  tested  the  effects  of  varying  the  mantle  Q  structure  and  find  that  the 
velocity  models  resulting  from  the  waveform  inversion  are  relatively  insensitive  to  reasonable 
variations  in  the  Q  structure.  Gomberg  Masters  (1988)  found  a  similar  result. 

We  next  tried  to  assess  the  effects  of  errors  in  our  data  and  assumptions  on  the  inversion 
results.  We  initially  assumed  the  CMT  source  mechanism  and  the  ISC  or  PDE  location  for  the 
events  in  inverting  the  regional  waveforms  but  found  that  in  several  cases  these  did  not  produce 
good  synthetic  fits  to  the  observed  seismograms.  Good  fits  could  be  obtained  after  slightly  altering 
the  CMT  source  mechanism  and/or  ISC  or  PDE  depth.  Both  the  CMT  source  mechanism  and 
the  final  mechanisms  and  depths  used  in  the  inversion  are  given  in  Table  2.  Errors  in  the  source 
mechanism  and  depth  will  map  into  the  earth  structure.  To  assess  the  effect  of  these  potential 
errors  we  have  performed  a  series  of  tests  in  which  we  held  all  but  one  of  the  source  parameters 
fixed  at  the  values  of  the  final  mechanism  in  Table  2,  varied  the  value  of  the  remaining  parameter 
by  a  reasonable  amount  and  reinverted  the  seismic  waveform  giving  a  new  earth  model.  We 
have  varied  the  strike,  slip,  and  dip  of  the  source  mechanism  by  ±10°  and  the  distance  and 
depth  by  ±10  km.  The  effects  of  depth  variations  are  so  pronounced  that  the  inversion  could 
not  converge  when  the  depth  varied  by  this  amount  and  we  conclude  from  this  that  the  regional 
surface  waveforms  provide  strong  constraints  on  the  focal  depth.  Of  the  remaining  parameters, 
the  largest  effects  are  for  variations  in  distance  and  dip.  The  bounds  shown  in  Figure  7c  represent 
the  largest  changes  in  the  inversion  earth  model  for  the  July  18,  1986  event,  resulting  from  these 
tests.  The  effect  of  these  variations  is  dependent  on  the  epicentral  distance  but  the  results  of 
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these  tests  suggest  that  for  the  magnitude  of  the  errors  in  distance,  strike,  dip,  and  slip  possible 
for  the  southern  African  events  we  have  used,  errors  in  velocities  could  be  as  large  as  ±0.1  km 
and  errors  in  density  could  be  as  large  as  ±0.05  g  cm~^. 

We  have  assumed  that  lateral  variations  in  earth  structure  along  the  propagation  path  are 
negligible  and  that  the  material  is  isotropic.  The  resulting  earth  model  could  be  biased  if  either 
of  these  assumptions  are  wrong.  Both  assumptions  are  difficult  to  test  with  the  data  available 
and  at  this  point  we  cannot  make  a  quantitative  statement  of  their  influence  on  our  result.  We 
have  shown  above  that  the  crust  of  southern  Africa  is  of  variable  thickness  but  in  general  is  41±6 
km  thick.  The  regional  surface  waveforms  are  well  fit  with  a  crust  of  almost  uniform  thickness  of 
41  km  over  southern  Africa.  The  major  variation  in  the  crustal  structure  results  from  variations 
in  the  seismic  velocities  in  the  shallow  (<  10  km)  crust.  The  nucleus  of  southern  Africa  consists 
of  the  Archaean  Zimbabwe  and  Kaapvaal  Cratons  which  are  dated  at  about  2.6  Ga,  and  the 
intervening  Limpopo  mobile  belt  of  similar  or  slightly  younger  age.  The  long  period  of  tectonic 
stability  of  the  region  might  suggest  that  lateral  variations  in  the  mantle  velocity  structure  of 
southern  Africa  are  also  relatively  minor.  However  Ballard  &  Pollack  (1987)  suggest  the  observed 
heat  flow  variations  in  southern  Africa  can  be  explained  by  a  thicker  lithospheric  root  beneath 
the  cratons  than  beneath  the  mobile  belts.  Some  global  models  show  thick  roots  under  the 
cratonic  regions  including  southern  Africa  but  there  is  no  consensus  on  this  issue.  A  somewhat 
thicker  lithospheric  root  beneath  the  cratons  is  also  suggested  by  the  distribution  of  diamond¬ 
bearing  kimberlites.  The  pressure  and  temperature  from  the  kimberlite  nodules  also  suggest  that 
the  lithospheric  thickness  varies  from  185  km  beneath  the  Kaapvaal  Craton  to  163  km  thick 
beneath  the  mobile  belts  (Finnerty  Boyd,  1987;  McKenzie,  1989).  Most  of  the  propagation 
paths  discussed  in  this  study  are  largely  confined  to  the  Archaean  core  but  they  do  extend  across 
portions  of  the  surrounding  mobile  belts  and  some  of  the  events  occur  in  the  southern  rift  (Fig. 
2).  While  lateral  variations  in  the  lithospheric  and  upper  mantle  structure  undoubtedly  occur,  the 
similarity  of  the  individual  path-average  earth  models  suggest  to  us  that  these  effects  are  small 
in  southern  Africa. 

Strong  anisotropy  in  the  lithosphere  or  upper  mantle  could  also  bias  our  results.  The  low 
frequency  surface  wave  observations  of  IMataf  et  al.  (1984)  show  the  upper  mantle  beneath 
southern  Africa  is  isotropic.  Vinnik  et  al.  (1995)  measured  shear  wave  splitting  from  SKS  and 
similar  phases  at  seven  sites  on  the  Kaapvaal  Craton  and  found  an  average  delay  of  0.9±0.5 
sec.  They  conclude  from  thermal  and  mineralogical  arguments  that  the  anisotropy  is  localized 
between  150  and  400  km  depth  requiring  2%  S-wave  anisotropy.  We  assume  that  the  lithosphere 
and  upper  mantle  are  isotropic  and  simultaneously  inverted  both  the  Love  and  Rayleigh  waves 
to  determine  the  path-average  earth  model.  For  the  frequency  range  of  our  data  we  found  no 
evidence  of  anisotropy  in  the  Love  and  Rayleigh  waves.  Thus,  while  we  cannot  make  quantitative 
statements  about  the  effects  of  lateral  variations  in  structure  or  anisotropy,  there  is  nothing 
apparent  in  our  data  to  indicate  that  these  effects  have  strongly  biased  our  inversion  results. 

We  have  next  determined  the  upper  mantle  compressional  wave  velocity  structure  by  forward 
modeling  the  Pn/  waveform  using  reflectivity  synthetics  (Fuchs  &  Muller,  1971).  We  first  con¬ 
verted  the  “blocky”  earth  model  determined  from  the  locked-mode  inversion  into  a  smoother 
earth  model  by  increasing  the  number  of  layers  and  changing  some  of  the  velocity  steps  in  the 
coarse  earth  model  into  gradients.  We  then  checked  to  see  that  these  smoothed  models  pro¬ 
duced  essentially  the  same  locked-mode  synthetic  seismograms  as  the  original  coarse  models. 
As  discussed  above,  we  used  generalized  ray  synthetic  seismograms  to  gain  some  insight  into 
the  rays  composing  the  initial  P-wave.  We  computed  reflectivity  synthetic  seismograms  for  the 
smoothed  earth  model  and  found  that  the  amplitude  of  the  synthetic  P-wave  first  arrival  was 
too  small  compared  to  the  observed  amplitude,  suggesting  that  there  was  insufficient  energy  in 
the  synthetic  mantle  turning  wave.  We  replaced  the  upper  mantle  P-wave  velocity  structure  of 
our  model  obtained  by  inversion  with  that  of  model  SACM06c  of  Clouser  &  Langston  (1990) 


and  found  that  this  produced  first  arrival  amplitudes  which  were  in  better  agreement  with  the 
observed  waveform  while  preserving  the  fit  of  the  surface  waves.  The  epicentral  distance 
of  the  July  18,  1986  event  to  SLR  is  similar  to  the  largest  epicentral  distances  of  events  in  the 
Clouser  Langston  (1990)  study.  When  we  modeled  events  at  larger  epicentral  distances  the 
synthetic  waveform  fits  were  not  as  good.  Through  a  trial-and-error  procedure  illustrated  in 
Figures  Bl-3,  we  finally  concluded  that  the  best  fit  over  the  entire  distance  range  of  our  data 
was  achieved  for  an  upper  mantle  compressional  wave  velocity  model  with  a  positive  gradient 
throughout  the  upper  mantle,  but  smaller  than  the  gradient  in  model  SACM06c. 

As  pointed  out  by  Clouser  Langston  (1990),  our  choice  of  upper  mantle  Qp  and  Qs  values 
affects  the  amplitudes  of  the  synthetic  waveforms.  If  we  had  chosen  higher  Q  values  then  a 
smaller  velocity  gradient  would  be  required  to  fit  the  P-wave  amplitude,  but  this  smaller  velocity 
gradient  would  result  in  a  degradation  of  the  travel  time  match.  It  does  not  seem  likely  that  a 
reasonable  change  in  the  Q  structure  would  remove  the  effects  of  a  LVZ  in  the  data. 

The  absence  of  a  P-wave  LVZ  leads  to  a  question  of  the  resolution  of  the  S-wave  LVZ  by  the 
regional'surface  wave  data.  We  performed  a  number  of  tests  to  ascertain  the  sensitivity  of  the 
regional  waveform  data  to  the  details  of  the  upper  mantle  structure.  The  Vs„  of  the  model  in 
Figure  B4a  is  constrained  by  the  Sn  observations  in  Figure  5  and  has  the  LVZ  replaced  by  a  positive 
gradient.  This  generates  synthetic  seismograms  which  have  larger  higher  mode  amplitudes  than 
those  seen  in  the  data,  suggesting  too  much  energy  in  the  synthetic  seismogram.  This  also  leads 
to  a  slight  mismatch  in  phase  resulting  in  the  higher  mode  and  the  low  frequency  fundamental 
mode  in  the  synthetic  leading  the  observed  phases.  We  found  it  impossible  to  find  S-wave  velocity 
models  which  fit  the  Sn  travel  times  from  Figure  5  but  which  did  not  contain  a  S-wave  LVZ. 

Stronger  evidence  for  the  existence  and  shape  of  the  S-wave  LVZ  comes  from  tests  on  more 
distant  events.  Figure  B4b  compares  observed  seismograms  for  the  July  24,  1991  earthquake 
(Event  10,  Table  1)  recorded  at  SUR  and  reflectivity  synthetics  computed  for  a  model  without  a 
S-wave  LVZ.  The  Sa  phase,  the  long  period  arrival  at  about  470  sec  and  primarily  seen  on  the 
vertical  component,  arises  from  the  coupling  of  energy  turning  in  the  upper  mantle  lid  and  in 
the  positive  gradient  at  the  base  of  the  LVZ  near  400  km  depth  (Priestley  et  al,  1980).  The  Sa 
phase  in  the  synthetic  seismogram  for  the  case  of  no  S-wave  LVZ  is  both  weaker  in  amplitude 
and  advanced  in  travel  time  compared  with  the  observed  Sa  phase.  On  the  other  hand.  Figure 
B4c  shows  that  for  a  somewhat  stronger  S-wave  LVZ,  the  Sa  phase  has  both  a  greater  amplitude 
and  delayed  travel  time  compared  to  the  observed  Sa  phase. 

Finally  we  tested  the  sensitivity  of  our  data  to  the  lid  thickness  by  comparing  the  observed 
seismograms  to  synthetic  seismograms  computed  for  models  with  the  thickness  of  the  lid  varying 
from  50  to  110  km  thickness.  Figure  B5a  shows  a  comparison  of  the  observed  waveforms  with 
synthetics  computed  for  an  earth  model  with  a  S-wave  LVZ  starting  at  90  km  depth.  The  higher 
mode  amplitude  is  less  than  that  of  the  observed  waveform,  suggesting  that  there  is  insufficient 
S-wave  energy  turning  from  the  upper  mantle  lid.  In  addition,  both  the  higher  mode  and  the 
low  frequency  fundamental  mode  are  delayed  with  respect  to  the  corresponding  observed  phases. 
Figure  B5b  shows  the  similar  comparison  for  a  lid  thickness  of  110  km,  that  is,  the  S-wave  LVZ 
starting  at  a  depth  of  150  km.  These  show  the  opposite  effect  -  larger  higher  mode  amplitudes 
and  advanced  travel  time  compared  to  the  observed  data.  Tests  for  smaller  changes  in  the  lid 
thickness  produced  changes  in  the  synthetic  waveforms  which  were  difficult  to  distinguish  from  the 
observed  waveform  and  we  therefore  feel  that  these  represent  maximum  bounds  on  the  thickness 
of  the  S-wave  lid. 

The  remaining  figures  of  this  appendix  (Figs.  B6-B16)  show  the  models  and  the  fit  of  the 
synthetic  and  observed  seismograms  for  the  remaining  eleven,  three-component  seismograms 
which  we  analyzed.  The  path  average  model  for  each  are  given  in  Tables  B1  to  Bll.  The  quality 
of  the  data  and  the  synthetic  -  observed  waveform  fits  are  discussed  in  the  captions  to  these 
figures.  83 


Table  B3  Earth  Model  for  March  9,  1989  -  SLR 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

kms  ^ 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

km 

P 

g  cm“^ 

0 

1 

4.51 

2.65 

2.29 

183 

17 

8.14 

4.39 

3.40 

1 

5 

6.38 

3.56 

2.76 

200 

17 

8.16 

4.36 

3.41 

6 

20 

6.48 

3.68 

2.94 

217 

17 

8.19 

4.34 

3.42 

26 

15 

6.73 

4.03 

3.08 

233 

17 

8.21 

4.31 

3.43 

41 

21 

8.00 

4.56 

3.28 

250 

27 

8.27 

4.35 

3.44 

63 

21 

8.01 

4.60 

3.30 

277 

27 

8.36 

4.45 

3.46 

84 

21 

8.01 

4.63 

3.32 

303 

27 

8.46 

4.55 

3.49 

105 

21 

8.02 

4.69 

3.35 

330 

27 

8.56 

4.65 

3.51 

126_ 

8 

8.02 

4.47 

3.37 

357 

27 

8.65 

4.75 

3.53 

134 

8 

8.04 

4.46 

3.37 

383 

27 

8.80 

4.90 

3.57 

142 

8 

8.07 

4.45 

3.38 

410 

31 

9.20 

5.10 

3.76 

150 

17 

8.09 

4.44 

3.38 

441 

19 

9.32 

5.17 

3.79 

167 

17 

8.11 

4.41 

0.00 

Table  B4  Earth  Model  for  March  10,  1989  -  SLR 


Depth 

km 

Thick. 

km 

Vp 

km  s“^ 

Vs 

kms“^ 

P 

g  cm~^ 

Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

km  s“^ 

P 

g  cm“^ 

0 

1 

4.51 

2.66 

2.30 

183 

17 

8.13 

4.39 

3.40 

1 

5 

6.34 

3.58 

2.77 

200 

17 

8.15 

4.36 

3.41 

6 

20 

6.42 

3.68 

2.95 

217 

17 

8.18 

4.34 

3.42 

26 

15 

6.71 

4.03 

3.07 

233 

17 

8.21 

4.31 

3.43 

41 

21 

7.99 

4.52 

3.27 

250 

27 

8.27 

4.35 

3.44 

63 

21 

7.99 

4.56 

3.29 

277 

27 

8.36 

4.45 

3.46 

84 

21 

8.00 

4.60 

3.31 

303 

27 

8.46 

4.55 

3.49 

105 

21 

8.02 

4.67 

3.35 

330 

27 

8.56 

4.65 

3.51 

126 

8 

8.02 

4.47 

3.37 

357 

27 

8.65 

4.75 

3.53 

134 

8 

8.03 

4.46 

3.37 

383 

27 

8.80 

4.90 

3.57 

142 

8 

8.05 

4.45 

3.38 

410 

31 

9.20 

5.10 

3.76 

150 

17 

8.07 

4.44 

3.38 

441 

19 

9.32 

5.17 

3.79 

167 

17 

8.10 

4.41 

0.00 
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Table  B5  Earth  Model  for  November  3,  1990  -  SLR 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

km 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

km 

P 

g  cm“^ 

0 

1 

4.49 

2.70 

2.09 

200 

20 

8.32 

4.43 

3.43 

1 

5 

5.18 

3.22 

2.43 

220 

20 

8.39 

4.41 

3.44 

6 

20 

6.35 

3.60 

3.11 

240 

13 

8.43 

4.41 

3.45 

26 

15 

6.65 

4.00 

3.18 

253 

13 

8.48 

4.43 

3.46 

41 

16 

8.05 

4.65 

3.35 

267 

13 

8.53 

4.45 

3.46 

57 

16 

8.07 

4.66 

3.35 

280 

22 

8.59 

4.49 

3.47 

73 

16 

8.09 

4.67 

3.36 

302 

22 

8.67 

4.56 

3.49 

89 

16 

8.10 

4.68 

3.37 

323 

22 

8.75 

4.63 

3.50 

105. 

16 

8.14 

4.69 

3.38 

345 

22 

8.83 

4.70 

3.51 

121 

20 

8.15 

4.50 

3.38 

367 

22 

8.91 

4.77 

3.53 

141 

20 

8.19 

4.48 

3.40 

388 

22 

9.03 

4.87 

3.55 

161 

20 

8.23 

4.47 

3.41 

410 

42 

9.36 

5.07 

3.76 

181 

20 

8.28 

4.45 

3.42 

452 

8 

9.50 

5.16 

3.81 

Table  B6  Earth  Model  for  July  24,  1991  -  SLR 


Depth 

km 

Thick. 

km 

Vp 

km  s“^ 

Vs 

kms~^ 

P 

g  cm~^ 

Depth 

km 

Thick. 

km 

Vp 

kms"^ 

Vs 

km 

P 

g  cm“^ 

0 

1 

5.00 

3.05 

2.50 

200 

20 

8.34 

4.43 

3.43 

1 

5 

6.20 

3.75 

3.00 

220 

20 

8.39 

4.41 

3.44 

6 

20 

6.20 

3.75 

3.00 

240 

13 

8.43 

4.41 

3.45 

26 

15 

6.75 

4.16 

3.16 

253 

13 

8.48 

4.43 

3.46 

41 

16 

8.11 

4.61 

3.35 

267 

13 

8.53 

4.45 

3.46 

57 

16 

8.13 

4.63 

3.35 

280 

22 

8.59 

4.49 

3.47 

73 

16 

8.15 

4.65 

3.36 

302 

22 

8.67 

4.56 

3.49 

89 

16 

8.17 

4.66 

3.36 

323 

22 

8.75 

4.63 

3.50 

105 

16 

8.21 

4.70 

3.37 

345 

22 

8.83 

4.70 

3.51 

121 

20 

8.22 

4.50 

3.38 

367 

22 

8.91 

4.77 

3.53 

141 

20 

8.25 

4.48 

3.39 

388 

22 

9.03 

4.87 

3.55 

161 

20 

8.28 

4.47 

3.40 

410 

42 

9.36 

5.07 

3.76 

181 

20 

8.31 

4.45 

3.41 

452 

8 

9.50 

5.16 

3.81 

86 


Table  B7  Earth  Model  for  July  24,  1991  -  SUR 


Depth 

km 

Thick. 

km 

Vp 

km  s“^ 

Vs 

kms"^ 

P 

g  cm~^ 

Depth 

km 

Thick. 

km 

Vp 

km  s“^ 

Vs 

km  s“^ 

P 

g  cm~^ 

0 

6 

5.81 

3.47 

2.69 

200 

20 

8.33 

4.48 

3.40 

6 

20 

5.97 

3.70 

2.79 

220 

20 

8.41 

4.48 

3.42 

26 

15 

6.66 

4.11 

3.00 

240 

28 

8.51 

4.51 

3.44 

41 

21 

8.07 

4.65 

3.32 

268 

28 

.8.60 

4.57 

3.46 

62 

21 

8.09 

4.68 

3.33 

297 

28 

8.70 

4.62 

3.48 

84 

21 

8.12 

4.70 

3.34 

325 

28 

8.79 

4.68 

3.50 

105 

21 

8.16 

4.76 

3.35 

353 

28 

8.89 

4.73 

3.52 

126 

25 

8.19 

4.59 

3.36 

382 

28 

9.03 

4.81 

3.55 

151. 

25 

8.22 

4.55 

3.37 

410 

42 

9.36 

5.01 

3.76 

175 

25 

8.27 

4.50 

3.39 

452 

8 

9.50 

5.10 

3.81 

Table  B8  Earth  Model  for  August  18,  1994  -  BOSA 


Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

kms“^ 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

km  s"^ 

P 

g  cm”^ 

0 

1 

5.91 

3.19 

2.50 

200 

20 

8.35 

4.33 

3.40 

1 

5 

6.28 

3.68 

2.84 

220 

20 

8.40 

4.31 

3.42 

6 

20 

6.30 

3.68 

2.86 

240 

13 

8.44 

4.30 

3.43 

26 

15 

6.74 

4.01 

2.97 

253 

13 

8.48 

4.30 

3.44 

41 

16 

8.10 

4.64 

3.30 

267 

13 

8.53 

4.30 

3.45 

57 

16 

8.12 

4.67 

3.31 

280 

22 

8.59 

4.34 

3.46 

73 

16 

8.14 

4.69 

3.32 

302 

22 

8.68 

4.43 

3.48 

89 

16 

8.16 

4.72 

3.33 

323 

22 

8.77 

4.52 

3.49 

105 

16 

8.20 

4.78 

3.34 

345 

22 

8.86 

4.60 

3.51 

121 

20 

8.20 

4.40 

3.35 

367 

22 

8.95 

4.69 

3.53 

141 

20 

8.24 

4.38 

3.36 

388 

22 

9.08 

4.82 

3.55 

161 

20 

8.28 

4.37 

3.38 

410 

42 

9.40 

5.02 

3.76 

181 

20 

8.31 

4.35 

3.39 

452 

8 

9.54 

5.16 

3.81 
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Table  B9  Earth  Model  for  August  18,  1994  -  LBTB 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

km 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

kms“^ 

Vs 

km 

P 

g  cm“^ 

0 

1 

5.90 

3.19 

2.51 

200 

20 

8.27 

4.33 

3.40 

1 

5 

6.24 

3.68 

2.85 

220 

20 

8.29 

4.31 

3.42 

6 

20 

6.27 

3.66 

2.84 

240 

13 

8.31 

4.30 

3.43 

26 

15 

6.73 

4.05 

2.96 

253 

13 

8.32 

4.30 

3.44 

41 

16 

8.15 

4.61 

3.28 

267 

13 

8.34 

4.30 

3.45 

57 

16 

8.16 

4.64 

3.29 

280 

22 

8.37 

4.33 

3.46 

73 

16 

8.17 

4.68 

3.31 

302 

22 

8.41 

4.38 

3.48 

89 

16 

8.18 

4.71 

3.32 

323 

22 

8.45 

4.44 

3.49 

,105 

16 

8.20 

4.78 

3.34 

345 

22 

8.50 

4.49 

3.51 

121 

20 

8.20 

4.40 

3.35 

367 

22 

8.54 

4.55 

3.53 

141 

20 

8.22 

4.38 

3.36 

388 

22 

8.60 

4.63 

3.55 

161 

20 

8.24 

4.37 

3.38 

410 

42 

9.21 

4.83 

3.76 

181 

20 

8.25 

4.35 

3.39 

452 

8 

9.36 

4.93 

3.81 

Table  BIO  Earth  Model  for  August  18,  1994  -  LSZ 


Depth 

km 

Thick. 

km 

Vp 

km 

Vs 

kms“^ 

P 

g  cm“^ 

Depth 

km 

Thick. 

km 

Vp 

km  s~^ 

Vs 

km 

P 

g  cm“^ 

0 

1 

5.92 

3.15 

2.50 

200 

20 

8.25 

4.34 

3.40 

1 

5 

6.29 

3.70 

2.84 

220 

20 

8.29 

4.31 

3.42 

6 

20 

6.35 

3.66 

2.85 

240 

13 

8.32 

4.30 

3.43 

26 

15 

6.75 

3.96 

2.95 

253 

13 

8.36 

4.30 

3.44 

41 

16 

8.07 

4.60 

3.24 

267 

13 

8.40 

4.30 

3.45 

57 

16 

8.09 

4.65 

3.26 

280 

22 

8.46 

4.35 

3.47 

73 

16 

8.11 

4.70 

3.29 

302 

22 

8.56 

4.44 

3.48 

89 

16 

8.12 

4.75 

3.31 

323 

22 

8.65 

4.54 

3.50 

105 

16 

8.16 

4.84 

3.36 

345 

22 

8.75 

4.63 

3.51 

121 

20 

8.17 

4.41 

3.36 

367 

22 

8.84 

4.73 

3.52 

141 

20 

8.19 

4.39 

3.37 

388 

22 

8.98 

4.87 

3.55 

161 

20 

8.21 

4.37 

3.38 

410 

42 

9.38 

5.07 

3.76 

181 

20 

8.23 

4.36 

3.39 

452 

8 

9.52 

5.16 

3.81 
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Table  Bll  Earth  Model  for  August  18,  1994  -  SDR 


Depth 

km 

Thick. 

km 

Vp 

km  s“^ 

Vs 

kms“^ 

P 

gcm"^ 

Depth 

km 

Thick. 

km 

Vp 

km  s”^ 

Vs 

km 

P 

g  cm~^ 

0 

1 

5.88 

3.16 

2.51 

200 

20 

8.18 

4.42 

3.41 

1 

5 

6.15 

3.61 

2.86 

220 

20 

8.20 

4.38 

3.42 

6 

20 

6.16 

3.67 

2.80 

240 

13 

8.23 

4.37 

3.43 

26 

15 

6.85 

4.18 

2.88 

253 

13 

8.28 

4.37 

3.44 

41 

16 

8.10 

4.52 

3.23 

267 

13 

8.34 

4.37 

3.45 

57 

16 

8.11 

4.58 

3.26 

280 

22 

8.42 

4.41 

3.46 

73 

16 

8.11 

4.63 

3.29 

302 

22 

8.53 

4.50 

3.48 

89 

16 

8.12 

4.69 

3.31 

323 

22 

8.64 

4.58 

3.49 

105 

16 

8.13 

4.80 

3.37 

345 

22 

8.75 

4.66 

3.51 

12'l 

20 

8.13 

4.53 

3.37 

367 

22 

8.86 

4.75 

3.52 

141 

20 

8.14 

4.50 

3.38 

388 

22 

9.03 

4.88 

3.55 

161 

20 

8.16 

4.48 

3.39 

410 

42 

9.43 

5.07 

3.76 

181 

20 

8.17 

4.45 

3.40 

452 

8 

9.57 

5.18 

3.80 
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APPENDIX  B  -  FIGURE  CAPTIONS 


Figure  Bl.  The  right  plot  shows  the  observed,  vertical  component  P-wave  of  Event  2  recorded 
at  SLR  (A  =  1038  km)  (solid  line)  compared  with  synthetic  seismograms  computed  for  the 
final  earth  model  for  this  event  from  Table  4  (dashed  line)  and  with  synthetic  seismograms 
computed  for  this  model  but  with  increasing  depths  to  the  top  of  the  P-wave  LVZ  (dotted 
line).  The  number  in  the  upper  left  hand  corner  of  each  seismogram  plot  denotes  the  depth 
to  the  top  of  the  P-wave  LVZ  for  the  model  used  to  compute  that  synthetic.  The  left  plot 
shows  the  P-wave  structure  of  the  final  earth  model  (solid  line)  and  of  LVZ  models  tested 
(dotted  lines).  Vp„  for  the  LVZ  models  tested  is  the  same  in  each  case  but  the  mantle  portion 
of  the  models  is  offset  to  the  left  in  the  velocity  plot  to  provide  a  clearer  comparison  with  our 
final  P-wave  velocity  model. 

Figure  B2.  This  is  a  similar  test  of  the  P-wave  LVZ  to  that  in  Figure  Bl  but  for  the  March  9, 
1989  earthquake  (Event  4,  Table  1)  recorded  at  SLR  (A  =  1479  km).  This  figure  is  in  the 
same  format  as  Figure  Bl. 

Figure  B3.  This  is  a  similar  test  of  the  P-wave  LVZ  to  that  in  Figure  Bl  but  for  the  July  24, 
1991  earthquake  (Event  10,  Table  1)  recorded  at  SUR  (A  =  2106  km).  This  figure  is  in  the 
same  format  as  Figure  Bl. 

Figure  B4.  (a)  Comparison  of  the  observed  seismograms  of  the  July  18,  1986  earthquake  (Event 
2,  Table  1)  recorded  at  SLR  (A  =  1038  km)  compared  to  synthetic  seismograms  computed 
for  the  final  earth  model  for  this  event  (Table  4)  (dashed  line)  and  the  same  earth  model 
but  with  the  S-wave  LVZ  replaced  with  a  constant  gradient  of  the  S-wave  lid.  (b)  Similar 
comparison  to  that  in  (a)  but  for  the  July  24,  1991  earthquake  (Event  10,  Table  1)  recorded 
at  SUR  (A  =  2106  km.  (c)  Comparison  of  the  observed  seismograms  for  the  same  event  as 
in  (b)  but  for  a  more  pronounced  S-wave  LVZ. 

Figure  B5.  Comparison  of  the  observed  seismograms  of  the  July  18,  1986  earthquake  (Event  2, 
Table  1)  recorded  at  SLR  (A  =  1038  km)  compared  to  synthetic  seismograms  computed  for 
the  final  earth  model  for  this  event  (Table  4)  (dashed  line)  and  for  the  same  earth  model  but 
with  the  S-wave  LVZ  starting  at  (a)  90  km  and  (b)  150  km  depth  (dotted  lines). 

Figure  B6.  July  7,  1983  -  SLR  (Event  1,  Table  1):  This  event  was  located  in  the  southern  part 
of  the  western  rift  and  was  recorded  at  SLR  (A  =  2038  km).  The  path  to  SLR  crosses 
the  Damara  and  Irumide  Mobile  Belts,  the  Zimbabwe  Craton,  the  Limpopo  Mobile  Belt,  and 
the  Kaapvaal  Craton.  About  50%  of  the  propagation  path  is  across  the  Archaean  core.  The 
epicenter  is  almost  due  north  of  SLR  so  the  observed  seismograms  are  nearly  naturally  rotated 
(N-S  radial,  E-W  transverse)  so  the  vertical,  north,  and  east  components  are  plotted.  The 
top  panel  of  the  figure  plots  the  whole  seismogram  (observed  -  solid  line,  synthetic  -  dotted 
line)  to  show  the  overall  waveform  fit;  the  middle  panel  plots  the  early  part  of  the  seismogram 
enlarged  to  show  the  details  of  the  lower  amplitude  P-  and  S-waves;  the  bottom  panel  plots 
the  final  earth  model  for  this  path  from  Table  Bl.  All  of  the  remaining  figures  of  this  appendix 
have  the  same  format. 

Figure  B7.  April  16,  1988  -  SLR  (Event  3,  Table  1):  This  event  occurred  on  the  northern 
boundary  of  the  Damara  Mobile  Belt  almost  due  north  of  SLR,  thus  the  seismograms  are 
naturally  rotated  to  radial  (N-S)  and  transverse  (E-W)  so  the  vertical,  north,  and  east  com¬ 
ponents  are  plotted.  The  1722  km  long  path  crosses  the  Damara  and  Irumide  Mobile  belts, 
the  Zimbabwe  Craton,  the  Limpopo  Mobile  belt,  and  the  Kaapvaal  Craton.  About  50%  of 
the  propagation  path  is  across  the  Archaean  core.  The  data  quality  for  this  event  is  poor 
compared  to  the  other  seismograms  used  in  this  study.  The  traces  are  affected  by  low  ampli- 
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tude,  long  period  noise  and  by  spiking  and  the  N-S  component  has  been  high-pass  filtered 
with  a  three  pole  Butterworth  filter  with  a  corner  at  0.05  Hz.  The  N-S  synthetic  was  filtered 
in  an  identical  manner.  This  noise  makes  the  body  wave  portion  of  the  seismograms  for  this 
event  unusable.  The  S-wave  velocity  model  shown  here  from  Table  B2  is  a  smoothed  version 
of  the  locked  mode  inversion  model. 

Figure  B8.  March  9,  1989  -  SLR  (Event  4,  Table  1):  This  event  is  located  in  the  Zambezi  Mobile 
Belt  and  the  1479  km  long  path  crosses  the  Zimbabwe  Craton,  the  Limpopo  Mobile  belt,  and 
the  Kaapvaal  Craton.  About  70%  of  the  propagation  path  is  across  the  Archaean  core.  The 
vertical  seismogram  is  well  fit  but  there  is  some  long  period  noise,  especially  apparent  on  the 
N-S  component.  The  model  plotted  in  the  lower  panel  is  from  Table  B3. 

Figure  B9.  March  10,  1989  -  SLR  (Event  5,  Table  1):  This  event  is  co-located  with  the 
March  9,  1989  event  but  was  larger  and  the  largest  amplitude  surface  waves  are  clipped.  The 
horizontal  seismograms  have  been  rotated  into  radial  and  transverse.  The  model  plotted  in 
the  bwer  panel  is  from  Table  B4. 

Figure  BIO.  November  3,  1990  -  SLR  (Event  8,  Table  1):  This  event  occurred  near  the  northeast 
boundary  of  the  Kaapvaal  Craton.  The  path  length  of  705  km  is  the  shortest  in  this  study, 
and  was  almost  entirely  confined  to  the  Kaapvaal  Craton.  The  E— W  component  seismogram 
has  been  high-pass  filtered  with  a  three  pole  Butterworth  filter  with  a  corner  at  0.05  Hz.  The 
E-W  synthetic  was  filtered  in  an  identical  manner.  Because  of  this  filtering  the  horizontal 
seismograms  have  not  been  rotated  into  radial  and  transverse.  The  model  plotted  in  the  lower 
panel  is  from  Table  B5. 

Figure  Bll.  July  24,  1991  -  SLR  (Event  10,  Table  1):  This  event  is  located  in  the  Zambezi 
Mobile  Belt  1071  km  northeast  of  SLR.  The  path  crosses  the  southern  portion  of  the  Zambezi 
Mobile  Belt,  the  Limpopo  Mobile  belt,  and  the  Kaapvaal  Craton.  About  60%  of  the  propaga¬ 
tion  path  is  confined  to  the  Archaean  core.  There  is  some  noise  and  spiking  apparent  on  all 
components.  Because  of  this  noise,  the  seismograms  have  not  been  rotated  and  the  vertical, 
north,  and  east  components  are  plotted.  The  model  plotted  in  the  lower  panel  is  from  Table 
B6. 

Figure  B12.  July  24,  1991  —  SUR  (Event  10,  Table  1):  This  is  the  same  event  as  in  Figure  Bll 
except  recorded  at  station  SUR  (A  =  2106  km).  This  propagation  path  is  similar  to  that  to 
station  SLR  but  also  crosses  the  Namaqua  Mobile  Belt.  About  50%  of  the  propagation  path 
is  confined  to  the  Archaean  core.  There  is  some  long  period  noise  especially  apparent  on  the 
horizontal  components.  Both  the  observed  and  synthetic  seismograms  have  been  band-pass 
filtered  with  a  three  pole  Butterworth  filter  with  corners  at  0.01  and  0.05  Hz.  The  horizontal 
seismograms  have  been  rotated  into  radial  (positive  away  from  the  source)  and  transverse 
(positive  clockwise  from  source).  The  model  plotted  in  the  lower  panel  is  from  Table  B7. 

Figure  B13.  August  18,  1994  ~  ROSA  (Event  12,  Table  1):  This  event  occurred  north  of  the 
Zambia  Craton  and  crossed  the  length  of  the  stable  region  of  southern  Africa.  About  60%  of 
the  propagation  path  to  BOSA  (A  =  2447  km)  is  confined  to  the  Archaean  core.  Both  the 
observed  and  synthetic  seismograms  have  been  low-pass  filtered  with  a  three  pole  Butterworth 
filter  with  a  corner  at  0.05  Hz.  The  horizontal  seismograms  have  been  rotated  into  radial  and 
transverse.  The  model  plotted  in  the  lower  panel  is  from  Table  B8. 

Figure  B14.  August  18,  1994  ~  LBTB  (Event  12,  Table  1):  This  is  the  same  event  as  in  Figure 
B13  except  recorded  at  station  LBTB  (A  =  2062  km).  About  50%  of  the  propagation  path 
is  confined  to  the  Archaean  core.  Both  the  observed  and  synthetic  seismograms  have  been 
low-pass  filtered  with  a  three  pole  Butterworth  filter  with  a  corner  at  0.05  Hz.  The  horizontal 
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seismograms  have  been  rotated  into  radial  and  transverse.  The  model  plotted  in  the  lower 
panel  is  from  Table  B9. 

Figure  B15.  August  18,  1994  ~  LSZ  (Event  12,  Table  1):  This  is  the  same  event  as  in  Figure  B13 
except  recorded  at  station  LSZ  (A  =  955  km).  Station  LSZ  is  located  in  the  Irumide  Mobile 
Belt  and  the  earth  model  derived  from  these  data  is  more  representative  of  the  Zambia  Craton. 
Both  the  observed  and  synthetic  seismograms  have  been  low-pass  filtered  with  a  three  pole 
Butterworth  filter  with  a  corner  at  0.05  Hz.  The  horizontal  seismograms  have  been  rotated 
into  radial  and  transverse.  The  model  plotted  in  the  lower  panel  is  from  Table  BIO. 

Figure  B16.  August  18,  1994  -  SUR  (Event  12,  Table  1):  This  is  the  same  event  as  in  Figure 
B13  except  recorded  at  station  SUR  (A  =  2996  km).  This  is  the  longest  path  length  of 
this  study  and  samples  the  length  of  the  stable  region  of  southern  Africa.  About  70%  of  the 
propagation  path  lies  within  the  Zambia,  Zimbabwe,  and  Kaapvaal  cratons  and  the  Limpopo 
Mobile  Belt.  Both  the  observed  and  synthetic  seismograms  have  been  low-pass  filtered  with 
a  three  pole  Butterworth  filter  with  a  corner  at  0.05  Hz.  The  horizontal  seismograms  have 
been  rotated  into  radial  and  transverse.  The  model  plotted  in  the  lower  panel  is  from  Table 
Bll. 
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FIGURE  B1 
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FIGURE  B5b 
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APPENDIX  C  Physical  Properties 

The  proportions  by  weight  of  the  end-members  in  the  various  phases  were  calculated  from 
the  molar  concentrations,  denoted  by  [  ],  of  the  elements  from  the  oxide  proportions  in  Table  Cl 
For  olivine  the  molar  concentrations  of  forsterite,  Fo,  and  fayalite,  Fa,  were  obtained  from 

Fo  = - _  Fa  - 

([Mg]  +  [Fe])’  "^-([Mg]  +  [Fe]) 

The  correspon^ng  expressions  for  pyroxenes  are  more  complicated.  If  [AI]>[Na],  the  proportions 
ot  jadeite,  Jd,  Ca-Tschermak,  CaTs,  and  Aegirine,  Ag,  are  first  obtained  from 

Jd  =  [Na],  CaTs  =  ([Al]  -  [Na])/2,  Ag  =  0 

and  the  molar  concentration  of  Ca  reduced  accordingly 


[Ca]  =  [Ca]  -  CaTs 


If  [Na]>[AI] 


Jd  =  [Al], 


and  the  molar  concentration  of  Fe  is 


Ag  =  [Na]  -  [Al], 


CaTs  =  0 


[Fe]  =  [Fe]-Ag 

The  iron  and  magnesium  are  distributed  between  Hd  &  Fs,  and  Di  &  En  by  defining 

p  _  |Mg] 

([Mg]  +  [Fe]) 

then 

Di  =  F[Ca],  En  =  ([Mg]-Di)/2 
Hd  =  (1  -  F’)[Ca],  Fs  =  ([Fe]  -  Hd)/2 

The  expressions  for  Feldspar  were  obtained  using 

An  =  [Ca],  Al  =  [Na],  Or  =  [K] 

Me'  in  accommodated  in  Sp  and 

Sp  =  [Al]/2,  Me  =  [Cr]/2 

Mf  =  [Mg]  -  (Sp  +  Me),  Mg  =  ([Fe]  -  Mf) /3 
He  =  Cr  =  0 

If,  however,  this  condition  is  not  satisfied,  but  [Mg]>[Cr]/2,  then 

Me  =  [Cr]/2,  Sp  =  [Mg]  -  Me 

He  =  [Al]/2-Sp,  Mg  =  ([Fe]  -  He)/3 
Mf  =  Cr  =  0 

Otherwise 

Me  =  [Mg],  He  =  [Al]/2 
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Cr  =  [Cr]/2  -  Me,  Mg  =  ([Fe]  -  He  -  Cr) /3 
Sp  =  Mf  =  0 

The  expressions  for  garnet  are  less  eomplieated.  If  F  is  defined  by 

f-  M 
(lAlJ  +  [Cr)) 

then 

Py  =  F[Mg],  A1  =  F[Fe],  Gr  =  F’[Ca] 

Kn  =  (l-F)[Mg],  Fe=(l-F’)[Fe],  Uv  =  (1  -  F)[Ca] 

Onee  the  molecular  proportions  of  the  end  members  had  been  obtained,  their  proportions  by 
weighTwere  calculated  using  the  molecular  weights  in  Table  Cl. 

The  densities  of  the  end  members  are  first  corrected  for  the  effect  of  temperature  using 

piT)  =  p  (l  -  ai(r  -  25)  -  a^iT  -  25f) 

where  02  =  0  except  for  olivine,  and  T  is  the  temperature  in  °C.  Then  the  density  and  seismic 
velocities  were  calculated  in  two  stages.  First  that  for  each  individual  phase  was  obtained  from 
the  proportions  of  the  end  members.  That  of  the  rock  as  a  whole  was  then  calculated  from  the 
proportions  of  each  phase.  The  expressions  are  the  same  for  each  calculation 

Up  =  EfilPi,  v''  =  -^v,uiPi 
2=1  2=1 


=  (y^  +  F^)/2 


where  /,  are  the  fractions  by  weight  of  the  end  members  or  of  the  mineral  phases,  and  Vi  is  either 
the  P  or  the  S  wave  velocity.  is  the  Voigt,  the  Reuss,  and  the  Voigt-Reuss-Hill 

average.  The  velocities  for  each  mineral  were  corrected  for  the  effect  of  temperature  and  pressure 


using 


F(P.7’)  =  V^(l  +  i(|I)(r-25)  +  i 


where  P  is  the  pressure  in  GPa.  When  the  final  P  and  S  wave  velocities  of  the  rock  had 
been  obtained,  the  rock  density  was  corrected  for  the  effect  of  pressure  using  the  bulk  modulus 
calculated  from  the  seismic  velocities 


p(T,P)  =  p(T,(l)  (1  +  P/(Vp^  - 

The  density  and  velocities  calculated  using  these  expressions  are  probably  not  accurate  to  better 
than  1-2%,  and  differences  in  density  and  velocity  caused  by  differences  in  mineralogy  to  better 
than  about  10%. 
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Table  Cl.  Mineral  Composition 


Si02 

Ti02 

AI2O3 

Cr203 

FeO 

MnO 

MgO 

CaO 

Na20 

K2O 

NiO 

P205 

olivine 

40.84 

0.02 

0.04 

0.04 

8.72 

0.50 

49.53 

0.06 

0.00 

0.00 

0.24 

0.00 

orthopyroxene 

57.32 

0.11 

0.89 

0.30 

5.35 

35.03 

0.12 

0.69 

0.17 

0.00 

0.01 

0.00 

clinopyroxene 

54.76 

0.21 

2.59 

1.63 

2.96 

0.10 

17.36 

18.34 

2.04 

0.01 

0.00 

0.00 

plagioclase 

44.79 

0.00 

34.81 

0.00 

0.27 

0.00 

0.03 

19.47 

0.64 

0.00 

0.00 

0.00 

A1  —  spinel 

0.00 

0.13 

65.55 

0.00 

11.94 

0.10 

22.28 

0.00 

0.00 

0.00 

0.00 

0.00 

Cr  —  spinel 

“0.12 

0.14 

10.28 

59.40 

17.06 

0.14 

12.62 

0.14 

0.00 

0.00 

0.00 

0.00 

garnet 

41.74 

0.43 

20.21 

3.99 

7.68 

0.37 

20.40 

5.11 

0.05 

0.00 

0.00 

0.00 

amphibole 

40.81 

5.25 

14.45 

0.03 

11.12 

0.13 

12.95 

10.69 

2.77 

1.81 

0.00 

0.00 

phlogopite 

42.79 

0.80 

15.25 

1.26 

3.55 

0.07 

25.90 

0.08 

0.70 

9.49 

0.13 

0.00 

apatite 

0.00 

0.00 

0.00 

0.00 

0.00 

0.07 

0.10 

55.84 

0.00 

0.00 

0.00 

42.05 

Calculated  modal  mineralogy  of  nodule  suites 


mg# 

ol 

opx 

cpx 

gar 

phi 

Cr  -  sp 

PKP 

93.1 

82.0 

7.77 

0.82 

0.36 

8.68 

0.19 

Depleted 

92.8 

78.2 

15.3 

1.67 

3.65 

0.64 

0.46 

GP  &  GPP 

91.4 

64.7 

23.3 

3.51 

6.70 

0.91 

0.75 

Fertile 

90.5 

62.5 

20.2 

7.08 

9.49 

0.41 

0.25 

Olivine 


V 


-6.4  X  10 


—5  1 


V 


dV 

dP 


y- 


13  X  10-^  GPa 


-1 


Forsterite 

Fayalite 

Mg2Si04 

Fe2Si04 

Fo 

Fa 

Mol.  wt. 

140.70 

203.78 

p  Mg  m“^ 

3.225 

4.392 

ai  °C-i 

2.84  X  10-® 

2.76  X  10- 

02  °C-2 

9.8  X  10-^ 

5.1  X  10"® 

Vp  km  s“^ 

8.57 

6.75 

Vs  km  s-i 

5.02 

3.41 
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Pyroxene 

V  (w)  =  V  (Ip)  = X 


Enstatite 

Ferrosilite 

Diopside 

Hedenbergite 

Jadeite 

Aegirine 

Ca  —  Tschermak 

Mg2Si206 

Fe2Si206 

CaMgSi206 

CaFeSi206 

NaAlSi206 

NaFeSi206 

CaAl2Si06 

En 

Fs 

Di 

Hd 

Jd 

Ag 

CaTs 

Mol.  wt. 

200.79 

263.86 

216.56 

248.10 

220.14 

284.97 

218.12 

p  Mg  m-3 

3.194 

3.998 

3.277 

3.632 

3.347 

3.576 

3.43 

a  °C-i 

2.9  X  10-3 

2.9  X  10-3 

3.2  X  10-5 

3.2  X  10-3 

3.0  X  10-3 

2.5  X  10-3 

2.6  X  10-3 

Vp  km  8“^ 

8.36 

6.90 

8.06 

7.22 

8.78 

7.32 

7.3 

Fs  km  s”^ 

.  -.4.99 

3.72 

4.77 

4.18 

4.82 

4.09 

4.1 

Feldspar 


1 

V  \dT 


=  -2.4  X  10 


°C-^  ^(^]=20x  10-3  GPa-^ 


V  \dP 


Anorthite 

Albite 

Orthoclase 

CaAl2Si208 

NaAlSiaOg 

KAlSisOg 

An 

A1 

Or 

Mol.  wt. 

278.2 

262.25 

278.25 

p  Mg  m-3 

2.762 

2.611 

2.56 

a  °C-i 

5.9  X  10-3 

5.7  X  10-3 

9.4  X  10-3 

Vp  km  s-^ 

7.29 

6.07 

5.91 

Fs  km  s-^ 

3.85 

3.40 

3.25 

Spinel 


V 


-20  X  10-3  °C-^  y  =  16  X  10-3  GPa-^ 


Spinel 

Magnesiocbromite 

Magnesioferrite 

Hercynite 

Chromite 

Magnetit' 

MgAl204 

MgCr204 

MgFe204 

FeAl204 

FeCr204 

FeFe204 

Sp 

Me 

Mf 

He 

Cr 

Mg 

Mol.  wt. 

142.27 

192.29 

199.99 

173.81 

223.84 

231.54 

p  Mg  m-3 

3.55 

4.43 

4.52 

4.40 

5.09 

5.20 

o  °C-^ 

2.7  X  10-3 

2.2  X  10-3 

2.2  X  10-3 

2.4  X  10-3 

1.8  X  10-3 

4.2  X  10- 

Vp  km 

9.25 

9.25 

9.25 

9.25 

9.25 

9.25 

Fs  km  s-^ 

5.01 

5.01 

5.01 

5.01 

5.01 

5.01 
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Garnet 


V  \dT) 


-4.8  X  10"®  °C"^ 


V 


=  8.5  X  10"®  GPa"^ 


Pyrope 

Almandine 

Grossular 

Knorringite 

— 

Mg3Al2Si30i2 

Fe3Al2Si30i2 

Ca3Al2Si30i2 

Mg3Cr2Si30i2 

Fe3Cr2Si30i2 

Py 

A1 

Gr 

Kn 

Fc 

Mol.  wt. 

403.15 

497.76 

450.45 

453.13 

547.78 

p  Mg  m~^ 

3.582 

4.318 

3.594 

3.852 

4.578 

2.6  X  10"® 

2.6  X  10"® 

2.3  X  10"® 

2.3  X  10"® 

2.0  X  10"® 

Vp  km  s~^ 

8.96 

8.42 

9.31 

8.50 

7.99 

V5  km.s"^ 

5.05 

4.68 

5.43 

4.79 

4.44 

Amphibole  Phlogopite  Apatite 


p  Mg  m  ^ 

3.05 

3.2 

3.2 

a  °C-i 

4  X  10"® 

4  X  10"® 

4  X  10' 

Vp  km  s"^ 

7.04 

5.55 

7.94 

V5  km  s~^ 

3.81 

2.88 

4.81 

V 


=  -9  X  10“®  °C"^ 


V 


30  X  10"®  GPa"^ 


r 


Uvarovite 

Ca3Cr2Si30i2 

Uv 

500.48 
3.848 
2.5  X  10"® 
8.88 
5.18 


117 


THOMAS  AHRENS 

SEISMOLOGICAL  LABORATORY  252-21 
CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
PASADENA,  CA  91 125 


RALPH  ALEWINE 
NTPO 

1901  N.  MOORE  STREET,  SUITE  609 
ARLINGTON,  VA  22209 


SHELTON  ALEXANDER 
PENNSYLVANIA  STATE  UNIVERSITY 
DEPARTMENT  OF  GEOSCIENCES 
537  DEIKE  BUILDING 
UNIVERSITY  PARK,  PA  16801 

RICHARD  BARDZELL 

ACIS 

DCI/ACIS 

WASHINGTON,  DC  20505 


DOUGLAS  BAUMGA'RDT 
ENSCO  INC. 

5400  PORT  ROYAL  ROAD 
SPRINGFIELD,  VA  22151 


WILLIAM  BENSON 

NAS/COS 

ROOM  HA372 

2001  WISCONSIN  AVE.  NW 

WASHINGTON,  DC  20007 

ROBERT  BLANDFORD 
AFTAC 

1300  N.  17TH  STREET 
SUITE  1450 

ARLINGTON,  VA  22209-2308 

RHETT  BUTLER 
IRIS 

1616  N.  FORT  MEYER  DRIVE 
SUITE  1050 

ARLINGTON,  VA  22209 

CATHERINE  DE  GROOT-HEDLIN 
SCRIPPS  INSTITUTION  OF  OCEANOGRAPHY 
UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 
INSTITUTE  OF  GEOPHYSICS  AND  PLANETARY  PHYSICS 
LA  JOLLA,  CA  92093 

SEAN  DORAN 

ACIS 

DCI/ACIS 

WASHINGTON ,  DC  20505 


RICHARD  J.  FANTEL 
BUREAU  OF  MINES 
DEPT  OF  INTERIOR  BLDG  20 
DENVER  FEDERAL  CENTER 
DENVER,  CO  80225 


MUAWIA  BARAZANGI 

INSTITUTE  FOR  THE  STUDY  OF  THE  CONTINENTS 
3126  SNEE  HALL 
CORNELL  UNIVERSITY 
ITHACA,  NY  14853 

T.G.  BARKER 

MAXWELL  TECHNOLOGIES 

P.O.  BOX  23558 

SAN  DIEGO,  CA  92123 


THERON  J.  BENNETT 
MAXWELL  TECHNOLOGIES 
1 1800  SUNRISE  VALLEY  DRIVE  SUITE  1212 
RESTON,  VA  22091 


JONATHAN  BERGER 
UNIVERSITY  OF  CA,  SAN  DIEGO 

SCRIPPS  INSTITUTION  OF  OCEANOGRAPHY  IGPP,  0225 
9500  GILMAN  DRIVE 
LA  JOLLA,  CA  92093-0225 

STEVEN  BRATT 
NTPO 

1901  N.  MOORE  STREET,  SUITE  609 
ARLINGTON,  VA  22209 


LESLIE  A.  CASEY 
DOE 

1000  INDEPENDENCE  AVE.  SW 
NN-40 

WASHINGTON,  DC  20585-0420 

STANLEY  DICKINSON 
AFOSR 

110  DUNCAN  AVENUE,  SUITE  B1 15 
BOLLING  AFB 

WASHINGTON,  D.C.  20332-001 
DIANE  I.  DOSER 

DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
THE  UNIVERSITY  OF  TEXAS  AT  EL  PASO 
EL  PASO,  TX  79968 


JOHN  FILSON 
ACIS/TMG/NTT 
ROOM  6T1 1  NHB 
WASHINGTON,  DC  20505 


1 


MARK  D.  FISK 

MISSION  RESEARCH  CORPORATION 
735  STATE  STREET 
P.O.  DRAWER  719 
SANTA  BARBARA,  C A  93 1 02-07 1 9 

LORI  GRANT 
MULTIMAX,  INC. 

31 1C  FOREST  AVE.  SUITE  3 
,PACIFIC  GROVE,  CA  93950 


ROBERT  GEIL 
DOE 

PALAIS  DES  NATIONS,  RM  D615 
GENEVA  10,  SWITZERLAND 


HENRY  GRAY 

SMU  STATISTICS  DEPARTMENT 
P.O.  BOX  750302 
DALLAS,  TX  75275-0302 


_  I.  N.  GUPTA 
’multimax,  INC. 

1441  MCCORMICK  DRIVE 
LARGO,  MD  20774 


JAMES  HAYES 
NSF 

4201  WILSON  BLVD.,  ROOM  785 
ARLINGTON,  VA  22230 


MICHAEL  HEDLIN 

UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 
SCRIPTS  INSTITUTION  OF  OCEANOGRAPHY  IGPP,  0225 
9500  GILMAN  DRIVE 
LA  JOLLA,  CA  92093-0225 

EUGENE  HERRIN 

SOUTHERN  METHODIST  UNIVERSITY 
DEPARTMENT  OF  GEOLOGICAL  SCIENCES 
DALLAS,  TX  75275-0395 


VINDELL  HSU 
HQ/AFTAC/TTR 
1030  S.  HIGHWAY  AlA 
PATRICK  AFB,  FL  32925-3002 


RONG-SONG  JIH 
PHILLIPS  LABORATORY 
EARTH  SCIENCES  DIVISION 
,  29  RANDOLPH  ROAD 

HANSCOM  AFB,  MA  0173 1-3010 

LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-200 
LIVERMORE,  CA  94551 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-221 
LIVERMORE,  CA  94551 


DAVID  HARKRIDER 
PHILLIPS  LABORATORY 
EARTH  SCIENCES  DIVISION 
29  RANDOLPH  ROAD 
HANSCOM  AFB,  MA  01731-3010 

THOMAS  HEARN 

NEW  MEXICO  STATE  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
LAS  CRUCES,  NM  88003 


DONALD  HELMBERGER 

CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 

DIVISION  OF  GEOLOGICAL  &  PLANETARY  SCIENCES 

SEISMOLOGICAL  LABORATORY 

PASADENA,  CA  91125 

ROBERT  HERRMANN 
ST.  LOUIS  UNIVERSITY 

DEPARTMENT  OF  EARTH  &  ATMOSPHERIC  SCIENCES 
3507  LACLEDE  AVENUE 
ST.  LOUIS,  MO  63103 

ANTHONY  lANNACCHIONE 
BUREAU  OF  MINES 
COCHRANE  MILL  ROAD 
POBOX  18070 

PITTSBURGH,  PA  15236-9986 
THOMAS  JORDAN 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 
EARTH,  ATMOSPHERIC  &  PLANETARY  SCIENCES 
77  MASSACHUSETTS  AVENUE,  54-918 
CAMBRIDGE,  MA  02139 

LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-207 
LIVERMORE,  CA  94551 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

LLNL 

POBOX  808,  MS  L- 175 
LIVERMORE,  CA  94551 


2 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN;  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-208 
LIVERMORE,  CA  94551 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-195 
LIVERMORE,  CA  9455 1 


THORNE  LAY 

UNIVERSITY  OF  CALIFORNIA,  SANTA  CRUZ 
EARTH  SCIENCES  DEPARTMENT 
EARTH  &  MARINE  SCIENCE  BUILDING 
SANTA  CRUZ,  CA  95064 

DONALD  A.  LINGER 
DNA 

6801  TELEGRAPH  ROAD 
ALEXANDRIA,  VA  22310 


LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
PO  BOX  1663,  MS  F665 
LOS  ALAMOS,  NM  87545 


LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
POBOX  1663,MSC335 
LOS  ALAMOS,  NM  87545 


KEITH  MCLAUGHLIN 
MAXWELL  TECHNOLOGIES 
P.O.  BOX  23558 
SAN  DIEGO,  CA  92123 


RICHARD  MORROW 

USACDA/IVI 

320  2 1ST  STREET,  N.W. 

WASHINGTON,  DC  20451 


JAMES  NI 

NEW  MEXICO  STATE  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
LAS  CRUCES,  NM  88003 


JOHN  ORCUTT 

INSTITUTE  OF  GEOPHYSICS  AND  PLANETARY  PHYSICS 
UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 
LA  JOLLA,  CA  92093 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-202 
LIVERMORE,  CA  94551 


LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 
ATTN;  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  808,  MS  L-205 
LIVERMORE,  CA  94551 


ANATOLI  L.  LEVSHIN 
DEPARTMENT  OF  PHYSICS 
UNIVERSITY  OF  COLORADO 
CAMPUS  BOX  390 
BOULDER,  CO  80309-0309 

LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
PO  BOX  1663,  MS  F659 
LOS  ALAMOS,  NM  87545 


LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
PO  BOX  1663,  MS  D460 
LOS  ALAMOS,  NM  87545 


GARY  MCCARTOR 

SOUTHERN  METHODIST  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
DALLAS,  TX  75275-0395 


BRIAN  MITCHELL 

DEPARTMENT  OF  EARTH  &  ATMOSPHERIC  SCIENCES 
ST.  LOUIS  UNIVERSITY 
3507  LACLEDE  AVENUE 
ST.  LOUIS,  MO  63103 

JOHN  MURPHY 
MAXWELL  TECHNOLOGIES 
1 1800  SUNRISE  VALLEY  DRIVE  SUITE  1212 
RESTON,  VA  22091 


CHARLES  ODDENINO 
BUREAU  OF  MINES 
810  7TH  ST.  NW 
WASHINGTON,  DC  20241 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K6-48 
RICHLAND,  WA  99352 


3 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K7-34 
RICHLAND,  WA  99352 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K6-40 
RICHLAND,  WA  99352 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K7-22 
^CHLAND,  WA  99352 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K6-84 
RICHLAND,  WA  99352 


FRANK  PILOTTE 
HQ/AFTAC/TT 
1030  S.  HIGHWAY  AIA 
PATRICK  AFB,  FL  32925-3002 


JAY  PULLI 

RADIX  SYSTEMS,  INC. 
6  TAFT  COURT 
ROCKVILLE,  MD  20850 


DAVID  RUSSELL 
HQ  AFTAC/TTR 
1030  SOUTH  HIGHWAY  AIA 
PATRICK  AFB,  FL  32925-3002 


SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5704 

MS  0979,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0979 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5791 

MS  0567,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0567 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5704 

MS  0655,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0655 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  6116 

MS  0750,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0750 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K5-72 
RICHLAND,  WA  99352 


PACIFIC  NORTHWEST  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 

PO  BOX  999,  MS  K5-12 
RICHLAND,  WA  99352 


KEITH  PRIESTLEY 
DEPARTMENT  OF  EARTH  SCIENCES 
UNIVERSITY  OF  CAMBRIDGE 
MADINGLEY  RISE,  MADINGLEY  ROAD 
CAMBRIDGE,  CB3  OEZ  UK 

PAUL  RICHARDS 
COLUMBIA  UNIVERSITY 
LAMONT-DOHERTY  EARTH  OBSERVATORY 
PALISADES,  NY  10964 


CHANDAN  SAIKIA 

WOOODWARD-CLYDE  FEDERAL  SERVICES 
566  EL  DORADO  ST.,  SUITE  100 
PASADENA,  CA  91 101-2560 


SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  41 15 

MS  0329,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0329 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  9311 

MS  1159,PO  BOX  5800 
ALBUQUERQUE,  NM  87185-1159 

SANDIA  NATIONAL  LABORATORY 
ATTN:  TECHNICAL  STAFF  (PLS  ROUTE) 
DEPT.  5736 

MS  0655,  PO  BOX  5800 
ALBUQUERQUE,  NM  87185-0655 

THOMAS  SERENO  JR. 

SCIENCE  APPLICATIONS  INTERNATIONAL 

CORPORATION 

10260  CAMPUS  POINT  DRIVE 

SAN  DIEGO,  CA  92121 


AVI  SHAPIRA 
SEISMOLOGY  DIVISION 

THE  INSTITUTE  FOR  PETROLEUM  RESEARCH  AND 
GEOPHYSICS 

P.O.B.  2286,  NOLON  58122  ISRAEL 

MATTHEW  SIBOL 
ENSCO,  INC. 

445  PINEDA  COURT 
MELBOURNE,  FL  32940 


JEFFRY  STEVENS 
MAXWELL  TECHNOLOGIES 
P.O.  BOX  23558 
SAN  DIEGO,  CA  92123 


DAVID  THOMAS 
ISEE 

29100  AURORA  ROAD 
CLEVELAND,  OH  44139 


LAWRENCE  TURNBULL 

ACIS 

DCI/ACIS 

WASHINGTON,  DC  20505 


FRANK  VERNON 

UNIVERSITY  OF  CALIFORNIA,  SAN  DIEGO 
SCRIPTS  INSTITUTION  OF  OCEANOGRAPHY  IGPP,  0225 
9500  GILMAN  DRIVE 
LA  JOLLA,  CA  92093-0225 

DANIEL  WEILL 
NSF 

EAR-785 

4201  WILSON  BLVD.,  ROOM  785 
ARLINGTON,  VA  22230 

RU  SHAN  WU 

UNIVERSITY  OF  CALIFORNIA  SANTA  CRUZ 
EARTH  SCIENCES  DEPT. 

1156  HIGH  STREET 
SANTA  CRUZ,  CA  95064 

JAMES  E.  ZOLLWEG 
BOISE  STATE  UNIVERSITY 
GEOSCIENCES  DEPT. 

1910  UNIVERSITY  DRIVE 
BOISE,  ID  83725 

DEFENSE  TECHNICAL  INFORMATION  CENTER 
8725  JOHN  J.  KINGMAN  ROAD 
FT  BELVOIR,  VA  22060-6218  (2  COPIES) 


ROBERT  SHUMWAY 
410MRAK  HALL 
DIVISION  OF  STATISTICS 
UNIVERSITY  OF  CALIFORNIA 
DAVIS,  CA  95616-8671 

DAVID  SIMPSON 
IRIS 

1 6 1 6  N.  FORT  MEYER  DRIVE 
SUITE  1050 

ARLINGTON,  VA  22209 

BRIAN  SULLIVAN 
BOSTON  COLLEGE 
INSITUTE  FOR  SPACE  RESEARCH 
140  COMMONWEALTH  AVENUE 
CHESTNUT  HILL,  MA  02167 

NAFI  TOKSOZ 

EARTH  RESOURCES  LABORATORY,  M.I.T. 
42  CARLTON  STREET,  E34-440 
CAMBRIDGE,  MA  02142 


GREG  VAN  DER  VINK 
IRIS 

1 6 1 6  N.  FORT  MEYER  DRIVE 
SUITE  1050 

ARLINGTON,  VA  22209 

TERRY  WALLACE 
UNIVERSITY  OF  ARIZONA 
DEPARTMENT  OF  GEOSCIENCES 
BUILDING  #77 
TUCSON,  AZ  85721 

JAMES  WHITCOMB 
NSF 

NSF/ISC  OPERATIONS/EAR-785 
4201  WILSON  BLVD.,  ROOM785 
ARLINGTON,  VA  22230 

JIAKANG  XIE 

COLUMBIA  UNIVERSITY 

LAMONT  DOHERTY  EARTH  OBSERVATORY 

ROUTE  9W  ' 

PALISADES,  NY  10964 

OFFICE  OF  THE  SECRETARY  OF  DEFENSE 
DDR&E 

WASHINGTON,  DC  20330 


TACTEC 

BATTELLE  MEMORIAL  INSTITUTE 
505  KING  AVENUE 

COLUMBUS,  OH  43201  (FINAL  REPORT) 


5 


PHILLIPS  LABORATORY 
ATTN:  XPG 
29  RANDOLPH  ROAD 
HANSCOM  AFB,  MA  01731-3010 


PHILLIPS  LABORATORY 
ATTN:  TSML 
5  WRIGHT  STREET 
SANSCOM  AFB,  MA  01731-3004 


PHILLIPS  LABORATORY 
ATTN:  GPE 
29  RANDOLPH  ROAD 
HANSCOM  AFB,  MA  01731-3010 


PHILLIPS  LABORATORY 
ATTN:  PL/SUL 
3550  ABERDEEN  AVE  SE 
KIRTLAND,  NM  871 17-5776  (2  COPIES) 


